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ABSTRACT 


A  hypothesis  is  made  that  the  Galerkin  Finite  Element 
Method  (GFEM)  offers  a  viable  option  to  the  traditional 
Finite  Difference  Method  (FDM)  for  numerical  weather  prediC' 
tion.   The  shallow  water  barotropic  primitive  equations  are 
the  forecast  equations  for  all  experiments.   The  hypothesis 
is  tested  by  observing  simple,  analytic,  atmospheric  wave 
propagation  on  uniform  and  variable  mesh  grids.   Second,  a 
strongly  forced  solution  simulating  small  scale  nonlinear 
interactions  is  evaluated  for  both  the  GFEM  and  FDM. 
Finally,  a  variable,  moving  grid  for  a  GFEM  model  is 
compared  to  a  uniform,  higher  resolution  GFEM  model  for  a 
strong  vortex  in  a  mean  flow.  The  GFEM  shows  a  better 
propagation  for  simple  atmospheric  waves  and  better  predic- 
tion to  a  forced  nonlinear  solution  than  the  FDM  model.   A 
moving  variable  grid  follows  an  area  of  strong  gradients 
while  not  generating  noise  in  the  transition  zone. 
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I.       INTRODUCTION 

Proper  simulation   of   atmospheric   flow   is   the   main 
objective   of   numerical    weather    prediction.      The    foundation 
of  numerical    weather   prediction   is    a    set    of   equations 
including  the  momentum,   thermodynamic,    moisture   and  conti- 
nuity  equations    for   modeling    the    atmosphere.      Through 
computer    simulation,    numerical    weather   prediction   predicts   a 
future    state   based    on    initial    conditions   describing    the 
atmosphere.      A   successful    forecast   model    must    include   small 
scale    processes.      Transports,    conversions    and    exchanges    of 
mass,    momentum    and    energy   occurring   on   the   small    scale 
represent   important   features   which   must   be   properly    simula- 
ted.     Feedbacks    from    the    proper    representation    of   these 
small    features   can    sometimes   markedly   influence   the   larger 
scale    solutions.      Representation    of   the    effects    of    small 
scale   processes   can    be    accomplished   directly   through 
increased    spatial    resolution.      Continued    increases    in    the 
spatial    resolution   will    eventually   allow    the   desired    process 
to    be    resolved    properly.      However,    a    doubling    of    resolution 
generally   requires    an    eight    fold    increase    in    computational 
effort.      The    value   of   increased    spatial    resolution    must 
ultimately   be   measured   by   its    contribution    to    the    overall 
forecast.      A   level    of   confidence    in    the    forecast    must    be 
achieved   that   the    process    is    resolved    near   that    particular 
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grid  resolution.   The  grid  resolution  could  be  uniformly 
fine  in  the  forecast  domain.  The  computational  effort  for 
uniform  fine  mesh  models  could  be  far  in  excess  of  available 
resources.   An  alternative  to  a  uniform  fine  mesh  is  a 
variable  mesh  where  the  fine  mesh  covers  only  regions  of 
interest  or  high  activity. 

The  conventional  numerical  weather  prediction  forecast 
scheme,  the  finite  difference  method,  approximates  the 
partial  differential  equations  with  a  truncated  Taylor 
series.   It  has  performed  admirably  when  forecasting  the 
larger  scales  of  motion.   Technological  improvements  in 
computing  power  coupled  with  better  understanding  of  the 
atmosphere  now  allow  the  smaller  scales  to  be  forecast. 

An  alternate  method  for  numerical  weather  prediction, 
the  GFEM  approximates  the  partial  differential  -equations 
while  minimizing  the  error  between  the  actual  equations  and 
their  approximation.   This  best  fit  logically  leads  one  to 
the  expectation  that  the  GFEM  will  better  model  the  smaller 
scales  than  the  finite  difference  schemes.   This  research 
will  demonstrate  practical  aspects  of  the  GFEM  rheoreti cal ly 
possible. 

In  this  dissertation,  the  GFEM  will  be  evaluated  to 
determine  its  potential  to  model  atmospheric  flow.   Equiva- 
lent GFEM  and  FDM  models  will  be  utilized  to  compare  the  two 
methods.   Simulations  of  small  scale  processes  explicitly 
resolved  on  uniform  and  variable  grids  will  be  investigated. 
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Rigorous  experiments  will  demonstrate  both  GFEM  and  FDM 
responses  for  forcing  near  the  grid  length  scale.   Finally, 
a  demonstration  will  be  made  of  a  GFEM  moving  variable  mesh 
which  moves  v^ith  the  small  scale  process  or  feature,  and 
thereby  allows  the  fine  mesh  to  resolve  the  highly  active 
region  initially  and  throughout  the  forecast  period. 
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II.   HYPOTHESIS 

Numerical  weather  prediction  has  steadily  improved  in 
the  simulation  of  atmospheric  flow.   Large  scale  flows  have 
been  adequately  represented  by  finite  difference  models  for 
several  decades.   Gal erki n-type  formulations  (Cullen,  1974b; 
Hinsman,  1975;  Staniforth  and  Mitchell,  1978;  Cullen  and 
Hall,  1979;  Staniforth  and  Daley,  1979;  MacPherson  and 
Aksel ,  1980;  and  Sasaki  and  Reddy,  1980)  have  been  shown  to 
be  competitive  with  finite  difference  models,  but  have  not 
shown  a  marked  improvement.  Comparing  the  current  opera- 
tional models  at  selected  large  computer  centers,  one  finds 
that  two  are  Galerkin  and  three  are  finite  difference 
(Galerkin:   National  Meteorological  Center  and  Canadian 
Meteorological  Center;  and  finite  difference:  Fleet 
Numerical  Oceanography  Center,  Air  Force  Global  Weather 
Center  and  European  Center  for  Medium  Range  Weather 
Forecasting.   However,  all  centers  have  ongoing  research 
with  Galerkin  models  and  there  are  indications  that  these 
will  give  better  long  range  forecasts.   The  dichotomy  arises 
because  the  Galerkin  applications  have  not  vindicated  them- 
selves with  a  marked  increase  in  accuracy,  but  rather  have 
shown  equivalent  accuracies.   Staniforth  and  Mitchel  (1977) 
stated  that  ultimately  the  best  gl  obal /hemi spheri c  models 
will  be  a  spectral  model.  The  spectral  model  is  based  on 
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the   Galerkin    procedure    and    the    use   of  trigonometric    func- 
tions  is   particularly    appealing    for    hemispheric    or   global 
grids.      However,    where    non-uniform    grids   are    required,    the 
GFEM    is    a    more    logical    choice. 

Marked    improvements    in    regional    forecasts   will    not   come 
from   better   large-scale   models,    but    rather   from    models    that 
can    simulate    smaller   atmospheric    features.      The   National 
Weather  Service  Limited   Fine   Mesh    Model    is   an    improvement 
over    the    hemispheric    model,    partly   because    it    has    higher 
resolution   and   can   resolve   smaller   phenomena.      These    smaller 
features    can    also    affect    the    large    scale    flow.      Numerical 
meteorologists   have   realized    this    for   years    and    have 
attempted    to    model    these    features. 

The  GFEM   hasthe    potenti.il    to    increase    efficiently    the 
spatial    resolution    for    the    purpose    of   simulating    accurately 
the    small-scale    processes.      '^  f    a    variable    mesh    is   to    be 
employed,    then    it    should   be    evaluated   by    proper   simulation 
of    a    small-scale    feature.       Ii    previous    research    (Staniforth 
and    Mitchell,    1978),    the    refined    grid    was    tested    by    pro- 
pagating   synoptic-scale   waves    into    the    finer    grid    and 
demonstrating    that    the    wave   could   move   into   the    finer   grid 
without    generation    of    significant    noise.      These   conditions 
represent    a    prerequisite.       If    synoptic-scale    waves    cannot 
move    freely    into    and    through    the    variable    grid,    then    advec- 
tion    interactions    will    not    occur    properly    in    the    fine    grid. 
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However,  one  must  also  show  that  a  smaller  scale  atmospheric 
feature  can  be  properly  resolved  or  developed  in  the  fine 
mesh  area.   This  will  be  a  milestone.   The  efforts  of 
Staniforth  and  Mitchell  (1978)  and  Staniforth  and  Daley 
(1979)  have  stressed  the  movement  of  synoptic  scale  features 
into  the  fine  mesh  area  but  have  not  demonstrated  improved 
resolution  of  smaller  scale  atmospheric  features  in  the  fine 
mesh  area. 

The  hypothesis  is  that  the  GFEM  is  a  viable  option  for 
numerical  weather  prediction  when  simulating  atmospheric 
flow  on  variable  grids.   Three  separate  features  of  the  GFEM 
will  be  explored.   Each  feature  will  establish  the  creden- 
tials of  the  GFEM  as  a  viable  option.   Each  feature  is 
intimately  related  to  the  proper  representation  of  a  small- 
scale  phenomenon.   The  cost  effectiveness  of  the  GFEM,  when 
evaluated  in  these  contexts,  will  provide  a  measure  of  the 
potential  contribution  of  the  method. 

A.    FEATURES 

The  following  three  specific  features  will  be  explored 
to  support  the  hypothesis. 

1.    Variable  Grid 

Investigations  of  a  suitable  alternative  to  the 
finite  difference  models  for  a  variable  grid  will  be 
performed.   Two  basic  subdivisions  are  avai 1 abl e--tri angul ar 
or  rectangular.   Staniforth  and  Mitchell  (1978)  employed  the 
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variable  rectangular  grid  as  shown  in  Figure  1.   The  grid 
contained  some  areas  of  finer  resolution  which  were  not  in 
the  verification  area.   Two  points  associated  with  having 
variable  resolution  in  peripheral  regions  arise: 

(a)  Unnecessary  computational  overhead  is  required; 
and 

(b)  Undesirable  phase  changes  occur  as  the  wave 
propagates  in  the  peripheral  regions. 

Older  (1981)  and  Woodward  (1981)  developed  c  transformation 
procedure  to  vary  smoothly  the  resolution  fcir  a  channel 
domain  from  a  coarse  to  a  fine  area  in  a  triangular  subdivi- 
sion. A  uniform  equilateral  subdivision  is  shown  in  Figure 
2.   The  use  of  triangles  allows  the  increased  resolution  to 
occur  only  where  desired  and  not  in  peripheral  regions. 
Two  possible  choices  of  triangles  include  the  right  and  the 
equilateral.   Hinsman  (1975)  utilized  equilc.teral  triangles 
while  Cullen  (1974b)  utilized  near-equilateral  triangles.   Both 
reported  excellent  wave  propagation.   Kelley  and  Williams 
(1976)  utilized  right  triangles  and  experienced  very  noisy 
solutions.   Woodward  (1981)  duplicated  Kelley  and  Williams 
effort  with  equilateral  triangles  and  found  a  major 
reduction  in  the  noise. 

The  differing  geometries  and  results  thus  far  reported 
mandate  a  further  review  of  triangular  subdivisions.  It  is 
not  obvious  which  subdivision  is  most  suitable.  A  distinct 
advantage  of  the  rectangular  subdivision  is  that  it  allows 
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Fig.  1.    Rectangular  subdivision  of  the  Northern  Hemisphere 
on  a  polar  stereographi c  projection  (Staniforth 
and  Mitchell,  1978). 
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X 


Fig.  2.   Triangular  subdivision  for  a  channel 
coordinates  (Woodward,  1981). 


in  Cartesian 
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algorithms  to  be  developed  which  take  full  advantage  of 
vector  processors.   However,  both  subdivisions  afford  the 
luxury  of  obtaining  a  variable  grid.   Therefore,  similar 
experiments  performed  on  each  type  subdivision  will  point 
out  the  advantages  and  disadvantages  of  each. 

2.    Numerical  Simulation  of  Physical  Processes 

The  variable  grid  chosen  as  the  most  suitable 
alternative  will  be  tested  while  simulating  a  physical 
process  to  show  that  a  small-scale  atmospheric  feature  can 
be  properly  portrayed  in  the  fine  mesh  area.   As  stated  in 
Chapter  II,  the  ability  of  variable  grids  to  resolve  small 
scale  atmospheric  features  has  not  been  rigorously  tested. 
The  ability  to  move  synoptic  scale  features  into  a  finer 
mesh  is  a  prerequisite  and  must  be  shown.  The  main  crux  of 
the  problem  is  what  is  happening  in  the  fine  mesh  area.  A 
proper  scheme  must  resolve  a  small  scale  feature  near  the 
smallest  grid  length.   Schoenstadt  (1980)  and  Williams 
(1981)  indicated  that  the  most  responsive  schemes  were  either 
a  staggered  finite  element  grid  with  primitive  variables  or 
an  unstaggered  finite  element  grid  with  vorti ci ty/di vergence 
formulation.   This  research  will  utilize  the  latter. 

The  simulated  process  will  be  that  of  a  mass  source 
analogous  to  that  found  in  the  upper  atmosphere  above  a 
hurricane  or  on  the  leading  side  of  a  strong  trough.   The 
source  will  appear  in  the  continuity  equation.   A  known  wave 
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form  will  be  assumed  and  the  required  source  to  support  the 
wave  form  will  be  analytically  derived.   The  expression  for 
the  wave  form  will  allow  control  of  the  scale  of  the  process 
so  that  a  near-grid  scale  phenomena  can  be  simulated. 
3.    Moving  Grid 

The  ability  to  move  the  fine  grid  si)  that  it 
remains  centered  on  an  atmospheric  circulation  will  also  be 
demonstrated.  This  capability  will  further  enhance  the 
applicability  of  the  GFEM  for  atmospheric  simulation.   If  it 
is  shown  that  small-scale  forcings  are    properly  reflected  in 
the  flow,  and  the  grid  can  be  moved  with  the  forcing 
disturbance,  then  the  GFEM  has  great  potential  for 
atmospheric  prediction. 
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III.   THE  DEVELOPMENT  OF  THE  GFEM 

The  hypothesis  in  Chapter  II  proclaims  the  GFEM  as  a 
viable  option  for  simulating  atmospheric  flow.   An  under- 
standing of  the  history  and  principles  of  the  GFEM  will  give 
further  insight  for  the  expected  improved  performance  as 
compared  to  the  conventional  FDM. 

The  origin  of  the  GFEM  can  be  traced  to  the  seventeenth 
century  work  by  Leonhard  Euler.   His  work  established  the 
branch  of  mathematics  known  as  calculus  of  variation.   He 
recognized  tiat  there  existed  a  partial  differential  equa- 
tion (PDE)  associated  with  the  minimization  of  a  functional. 
This  POE  has  been  appropriately  named  the  Eul  er-Lagrange 
equation.   Solving  the  PDE  was  equivalent  to  determining  a 
stationary  value  of  the  functional.   In  the  late  nineteenth 
century.  Lord  Rayleigh  furthered  the  variational  calculus  by 
representing  the  dependent  variable  as  a  mathematical 
expression,  typically  as  a  power  series.   The  stationarity 
of  the  solution  allowed  determination  of  the  unknown  coeffi- 
cients of  the  power  series.   His  work  was  generalized  in  the 
early  twentieth  century  by  W.  Ritz  and  the  procedure  has 
since  been  referred  to  as  the  Rayl ei gh-Ri t z  method.   A 
limitation  in  the  solution  by  the  Rayleigh-Ritz  method  is 
the  fullness  of  the  matrix  to  be  inverted.   Shortly  after 
Ritz  completed  his  work,  Galerkin  developed  a  procedure 
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using  weighted  residuals  which  had  a  wider  application  than 
classical  variational  calculus.   The  residual  is  made  ortho- 
gonal to  a  function,  called  the  test  function,  and  thereby  is 
minimized.   While  the  Rayl ei gh-Ri tz  procedure  applies  to  the 
functional  of  an  associated  Eul er-Lagrange  equation,  the 
Galerkin  procedure  pertains  to  any  PDE.   The  methods  produce 
identical  results  when  applied  to  equivalent  extremum 
probl ems . 

The  Galerkin  procedure  was  unknowingly  well  established 
by  the  end  of  World  War  II.   The  rapid  technological 
advances  made  during  and  after  the  war,  coupled  with  the 
emergence  of  the  computer,  led  the  aircraft  industry  to 
develop  new  numerical  procedures  for  solving  stress  problems 
in  aircraft  design.   A  successful  technique  was  developed 
and  mathematicians  realized,  long  after  the  fact,  that  the 
Galerkin  procedure  was  being  utilized. 

The  GFEM  is  a  commonly  used  subset  of  the  set  of 
Galerkin  procedures.   After  subdivision  into  a  set  of  ele- 
ments, the  domain  resembles  a  completed  jigsaw  puzzle,  and 
hence  leads  to  the  terminology  "finite  elements."   Figures  1 
and  2  are  samples  of  such  subdivisions.   The  dependent 
variables  of  the  PDE  are  represented  as  a  linear  combination 
of  known  functions,  which  are  usually  low  order  polynomials. 
The  same  function  is  employed  as  the  test  function.   When 
substituted  into  the  PDE,  these  approximations  leave  a 
residual.   Minimizing  the  residual  completes  the  procedure. 
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The  functions  are  globally  zero  except  where  the  dependent 
variables  are  defined  near  each  individual  node.   The  pro- 
duct of  functions  remains  zero  except  where  the  representing 
and  test  functions  are  both  non-zero.   This  makes  the  method 
attractive  for  computer  implementation.   It  removes  the 
matrix  "fullness"  problem  found  with  the  Rayl ei gh-Ri tz 
approach.   The  minimization  process  lies  at  the  heart  of  the 
expected  improved  performance.   Each  term  of  the  POE  has 
been  simultaneously   approximated  and  the  error  in  those 
approximations  has  been  minimized. 
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IV.   MODEL  DESCRIPTION 

Two  different  barotropic  shallow-water  models  will  be 
employed  to  test  the  hypothesis  in  Chapter  II.   Each  model 
will  have  the  identical  domain,  boundary  and  initial 
conditions.   The  difference  will  lay  with  the  partial 
differential  equation  approximation.   In  one  model,  a 
Galerkin  approximation  to  the  partial  differential  equation 
will  be  used,  while  the  other  model  will  use  a  finite 
difference  approximation.   Two  versions  of  the  Galerkin 
model  will  be  evaluated.   The  first  version  will  have  trian- 
gular subdivisions  and  basis  functions,  while  the  other  will 
have  rectangular  subdivisions  and  basis  functions. 

A.    GALERKIN  FINITE  ELEMENT  MODEL 

The  system  of  equations  referred  to  as  the  shallow- 
water  equations  consists  of  three  equations  with  three 
forecast  variables  <p  y    u  and  v.   The  equations  are  written 


3t    3x    ay   ^^3x   ay^ 


9u    lu    lu  _  ^    ii  ^  Q 
at    ax    ay       ax 


(4-2) 


IV    IV    av  ^  ^    M  =  0 
at    ax    ay       ay 


(4-3) 
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Here  <j)  is  the  geopotential  height, 

u  is  the  east/west  component  of  the  wind, 

V  is  the  north/south  component  of  the  wind,  and 

f  is  the  Coriolis  parameter 

By  expanding  cj)  into  a  mean  ($)  and  a  deviation  (cj)')  the 

equations  can  be  written 


(4-4) 


du    _^  _3i.'  ,  9K 


3t 


3X 


ax 


Q  =  0 


(4-5) 


at   ay    ay    ^ 


(4-6) 


'y    3y 

here   D  is  the  divergence 

K  is  the  kinetic  energy  (per  unit  mass),  and 

Q  is  the  absolute  vorticity. 
The  primes  will  be  dropped  for  the  rest  of  the  paptr  for 
cl  ari  ty . 

Cullen  and  Hall  (1979)  showed  that  the  accuracy  of  the 
GFEM  solution  was  better  for  the  vorti ci ty-di vergence 
formulation  of  the  shallow-water  equations  than  for  an 
increase  in  resolution  with  the  primitive  formulation. 
Williams  and  Schoenstadt  (1980)  noted  that  staggered 
variable  formulation  of  the  primitive  equations  and  the 
unstaggered  vorti ci ty-di vergence  formulation  gave  the  best 
treatment  of  geostrophic  adjustment  for  small-scale 
features. 
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The   vorti ci ty-d1 vergence   form   of  the   shallow-water 
equations   also   allows   the   use   of   a   semi -i mpl i ci t   time 
scheme.      This    scheme   artificially   slows   the   propagation 
speed   of  the   fastest   gravity   waves,    which   allows   a   much 
larger   time    step    than    one    could    expect    for   a   normal    Courant- 
Fredri ch-Lewy    (CFL)    stability   criterion.      This    scheme    thus 
offsets    some   of   the    extra   computational    expense    required    to 
solve   the   system    of   equations    assembled    at    each   time    step. 

The    vorti ci ty/di vergence    form    of   the   equations    is 

II  +    4)0    +   ^(u6)    +  ^;:(vcD)    =    0  (4-7) 


^  +   -^(uQ)    +  T-^-(vQ)    =    0 


3t 


3X 


yy 


(4-8) 


1^  +    v%    +    V^K    -   ^(vQ)    +   ^(uQ)    =    0  (4-9) 

Here  v^  is  the  Laplacian  operctor  and  c  is  the  {—    -tSt) 

0  X    o  y 

relative  vorticity. 

The  velocity  can  be  writ':en  a.i   the  sum  of  the  rotational  and 

irrotational  components  as 


V  =  V,  +  V 


where  V,  =  KXV'^  and  V   =  Vx 


The  equations  can  be  rewritten  using 

'^    3X    3y 


and 


,  =  IV  .  iu  ,  ^2 
^   ax   3y 
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resulting  in 


at 


J^M    -   ^  (V*) 


(4-10) 


9  o2, 
3t ^  ^  = 


^(uQ)  -  ^(vQ) 


(4-11) 


(4-12) 


|-rv^x  +  v^(|,  =  -^CvQ)  -  — (uO)  -  -^-  ii  -  -L  3J< 
^t  3x^^^^   ay^"^^   3x  3x   ay  37 

Equation  (4-12)  can  be  further  manipulated 

The  domain  of  integration  is  a  channel  with  east-west 
periodicity.   The  boundary  condition  at  a  wall  is 

V  .  N  =  0 
where  N  is  an  outward  pointing  normal  vector.   Along  the 
northern  and  southern  walls,  the  v  component  is  equal  to 
zero,  so  that  the  v  equation  of  motion  (4-3)  reduces  to 

If '  -f" 

The  zonal  and  meridional  components  of  the  wind  can  be 
written  as  non-divergent  and  irrotational  components  as 


and 


u  =  .  11  +  iX 
ay   3x 


ax   ay 


Then,  along  the  north/south  walls  where  v  equals  zero,  the 
boundary  condition  is 


ax   ay 


(4-14) 
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The  above  condition  is  imposed  by  setting 

ij;  =  a  constant 
when  solving  the  vorticity  equation  and 


(4-15) 


IX  =  0 

3y 


(4-16) 


when  solving  the  divergence  equation.   This  is  an    overspe- 
cification  but  (4-14)  would  be  difficult  to  apply.   The 
initial  conditions  which  will  be  presented  later  are 
specifically  selected  to  satisfy  (4-15)  and  (4-16), 

The  semi -i mpl i ci t  scheme  is  implemented  by  evaluating 
all  the  terms  on  the  left  hand  side  of  the  equations  as  an 
average  at  time  levels  (t  +  At)  and  (t  -  At)  or  with  a 
centered  time  difference  as  appropriate.   All  the  terns  on 
the  right  hand  side  are  evaluated  at  time  level  t.   The 
equations  become 

||-  +  $(v^(t  +  At)  +  V^Ct-At))  =  -  ^^-(uQ)  -  ^(vQ)     (4-17) 


v2  11  = 
St 


^(uQ)  -  ^(vQ) 


V^(|f  +  (j)(t  +  At)  +  <(.(t-At))  =  -^((VQ) 


3X^ 


(4-18) 


-i_((L,Q)  +  iii) 

ay       dy' 


(4-19) 


The  divergence  equation  (4-19)  can  be  solved  for  x  at 


(t  +  At)and  substituted  into  the  equation  (4-17)  to  yield 


V  ^ 


$(At) 


3^((vQ)  -If)  -^((UQ)  .|f) 


$  (  A  t ) 


(4-20) 
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where  the  overbar  denotes  an  average  of  t  +  At  and  t  -  At. 
The  system  of  three  equations  in  three  unknowns  has  been 
reduced  to  two  Poisson  equations  and  one  Helmholtz  equation 
to  be  solved  at  each  time  step. 

The  solution  orocedure  involves  solving  the  <^   equation 
(4-20)  for  a  new  i.   The  divergence  equation  (4-19)  is  then 
solved  for   ^  "*■  Tt  ^"^  ^^    subtraction  for  y|-.   Finally,  the 
vorticity  equation  (4-18)  is  solved  for  v±-.   There  are  two 
options  as  to  whici  variables  will  be  history -carrying: 
either  <}),  ii   and  x  or  <!>,  u  and  v.   The  choice  was  made  to  use 
<j),  u  and  V  as  history-carrying  variables.   They  are  updated 
after  each  time  sttjp  following  Staniforth  and  Mitchell 
(1977)  by 

(J)(t  +  At)  =  2j  -  ({.(t-At)  (^-2T  ) 

u(t.At)  =  2At(^fX-  ^M)  .  u(t-At)        (4-22) 

v(t.:t)  -  2AtC^|X.^M)  .  v(t-At)        (4-23) 

Implementation  of  the  GFEM  is  accomplished  as  described 
in  Chapter  III.   The  north  and  south  latitudes  are  input 
parameters,  and  the  north/south  distance  is  subdivided  into 
N  equal  parts,  where  N  is  also  an  input  parameter.   The 
east/west  increment  is  a  function  of  the  north/south  incre- 
ment, as  explained  in  the  section  describing  the  individual 
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experiments.   The  Coriolis  parameter  will  be  a  constant 
equal  to  the  mid-channel  value. 

The  triangular  and  rectangular  domain  subdivisions  will 
be  demonstrated  in  the  experiment  sections.   An  appropriate 
approximating  function  for  the  rectangular  subdivision  is  a 
bilinear  function.   In  either  case,  the  forecast  equations 
in  Galerkin  form  are 

o)V,.  =  /(^((vQ),V,  -  :^ 
(At) 


/^^V^- w^^•=M^^^^^^^•■-^•^^^^^ 


-  /^^(^Q^j^j '  ^^J^^)^• '  /^(^^*^)j^j  ^  jj^^'hh^'i 


-I^r^^'-''^h'i  'jit'S^'-^'^'^h 


(4-24) 


fv^  11  .v.V.  =  -  f(— (uQ)  .V.)V.  -  f(:r-CvQ)  .V.)V. 


(4-25) 


/v^C^x^.v.  .J.V.IV.  =/C^CCvQ).V 


—  .V.  ))V. 
3x  J  j^^  i 


-/ 


Ct^(CuQ)  .v.  +  —  .V.))V. 


(4-26) 


Here  the  integral  sign  implies  an  area  integral  over  the 
domain,  the  j  subscript  denotes  Einstein  summation  for  the 
dependent  variables  and  the  i  subscript  is  the  ith  nodal 
equati  on . 
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The  strength  of  the  GFEM  is  that  the  spatial  deriva- 
tives of  the  dependent  variables  become  derivatives  of  the 
known  approximating  function  and,  therefore,  are  known 
exactly  at  each  node.  When  using  piecewise  continuous 
linear  functions,  the  first  derivatives  are  piecewise 
discontinuous  and  second  derivatives  are  not  defined. 
Therefore,  the  second  derivatives  must  be  handled  in  a 
different  fashion, 

2  .2 

3x-   -  ■    3y' 

2 

3X 


|v^.v.)=/(ii,(v.v^).-l!(v.v^)) 


/5XJ  •!    yax^  8x 


/"i__V.v  =  /"J-.-iiiV,-)  .  riV.  iV. 

J    3y2    J     i         J3y-3yJ     '^        j    3yJ     3y^ 


J   3y     ;  3y  J  3y  i 


where  the  <p  implies  a  line  integral  along  the  east/west 
or  north/south  boundaries.  The  east/west  line  integrals 
are  zero  because  of  peri odi ci ty, but  the  north/south  line 
integrals  add  additional  terms  to  the  equations  and  are 
satisfied  by  the  boundary  conditions.   The  Laplacian  terms 
become 
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/'S^  -M^■^^  -fjv'sPi  ^^i7^^ 


This    integration-by-parts    procedure    allows    the    use   of   a 
linear   function   while   approximating    second   derivatives.      The 
final    form   of   the    forecast    equations  is 


J     3X"^     3X 


J    ay"    3y 


$(At) 


/ 


9V 


(^((^'"j^  -  MiJ"^ 


/■,9*.    3V.    3V.     .     i*.     3V.     3V.^  f(l..n\       9V.,,, 

-  ;<3tJ  ^J  i7^  *  it  J  iyj  iy^ '  "  "  /''"'''j  i^J'^ 


-/ 


(Cvq).  I^j)  V. 


(ll^A    3V.    3V.    ^  ix-    3V .    aV.     .    -    w    J    N 

-     /(~rJ    — 1    — 1    +   ~rl    — 1    — 1    +    (j).V.V.)    = 

i^9t^    ax^    3X  3t^    3^^    3y  ^J    J    V 


/ 


^^^^^QjjVj  -  s- 17^■^^^• 


-/(t^   C(uQ)  .v.    +    K.   — j))V. 


(4-27) 


(4-28) 


(4-29) 
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Here  the  Helmholtz  equation  (4-27)  has  utilized 


and 


V^X  =  DIVERGENCE  =  ~  +  ~ 

a  A     ay 


ii  =  -  u  f 
9y       0 

al ong  the  wal 1 s. 

The  line  integral  along  the  north  and  south  boundary 
has  been  dropped  from  the  vorticity  ecuation  (4-28),  since 
the  value  of  ^    on  the  boundaries  is  a  constant  and  therefore 
the  value  of  -A:   is  zero.   The  line  integral  along  the  north 
and  south  boundary  has  been  dropped  from  the  divergence 
equation  (4-29),  because  the  value  of  the  normal  derivative 
along  the  north/south  walls  is  zero,  cs  stipulated  by  the 
boundary  conditions.   The  initial  conditions  will  also 
satisfy  the  condition  that  the  normal  derivative  of  x  is 
zero  along  the  north/south  walls. 

B.    FINITE  DIFFERENCE  MODEL 

The  comparison  model  to  the  GFEM  model  is  an  adaptation 
of  the  staggered,  primitive-equation  model  as  described  in 
Section  7-4  of  Haltiner  and  Williams  (1980).   The  equations 
are  in  the  flux  form  and  employ  Scheme  C  as  shown  in  Section 
7-3  of  Haltiner  and  Williams  (1980).   The  baroclinic  model 
described  there  has  been  coded  to  perform  in  a  barotropic 
mode. 
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Since  the  GFEM  is  written  in  differentiated  form, 
vorticity  and  divergence  have  been  diagnosed  from  the  primi- 
tive variables  in  the  finite  difference  model  to  allow 
appropriate  comparisons. 

C.    NUMERICAL  METHODS 

Although  no  attempt  is  made  to  optimize  computational 
efficiency  in  this  research  model,  the  GFEM  must  utilize 
efficient  numerical  techniques  to  be  considered  a  viable 
option  for  numerical  weather  prediction.   Various  solution 
procedures  are  available  for  the  GFEM  system  of  equations. 
Successive  over-relaxation  was  employed  during  the  early 
research  stages.   Several  problems  arose  which  indicated  a 
need  for  a  better  solution  procedure.   The  successive  over- 
relaxation  procedure  resulted  in  a  bias  if  the  sweeping 
during  each  pass  was  in  the  same  direction.   Alternating  the 
direction  alleviated  that  particular  problem,  but  the  number 
of  passes  required  to  achieve  a  desired  level  of  accuracy 
became  prohibitive.   A  second  approach  employed  a  direct 
solver  using  a  Gaussian  elimination  procedure.   The  matrices 
from  the  Galerkin  procedure  are  decomposed  into  upper  and 
lower  block  tri-diagonal  matrices.   A  preprocessing,  repre- 
senting the  forward  substitution  stage,  can  be  done  once. 
The  back  substitution  must  be  done  each  time  a  solution  is 
desired.   The  coefficients  necessary  for  this  back  substitu- 
tion are    stored  in  an  efficient  manner.   The  particular 
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algorithm  is  called  a  "skyline"  solver  and  is  described  in 
detail  in  Bathe  and  Wilson  (1976).   Skyline  refers  to  the 
compact  method  of  storing  only  those  coefficients  required. 
This  method  has  the  desired  level  of  accuracy  and  a  high 
degree  of  computational  efficiency.  Staniforth  and  Mitchell 
(1977)  employed  another  direct  solver  method  called,  the 
conjugate-  gradient  method.   This  method  is  yery   computa- 
tionally efficient,  even  when  the  non-constant  coefficient 
Helmholtz  equation  is  involved.   The  theoretical  oper'ation 
count  indicates  that  for  very  large  domains  with  many 
degrees  of  freedom,  the  conjugate-gradient  method  wil  be 
much  better  than  the  other  direct  solver  mentioned  a^ove. 
The  Gauss  elimination  procedure  was  utilized  instead  of  the 
conjugate-gradient  method  for  reasons  of  expediency. 

Numerical  integration  of  the  three  forecast  equations 
involves  solving  first  a  Helmholtz  equation  for  cb ,  s€;cond  a 
Poisson  problem  for  ip    and  finally  a  Poisson  problem  for  x« 
The  boundary  conditions  and  the  equations  are    well  posed  for 
the  first  two  equations.  However,  the  normal  derivative  of 
X  is  equal  to  zero  for  the  last  Poisson  equation  and 
requires  special  attention.   The  solution  plus  a  constant  is 
also  a  solution.   Since  only  derivatives  of  x  ^re  required, 
the  shape  of  the  field  is  much  more  important  than  the 
actual  value  of  x*   To  avoid  computer  round  off  errors,  the 
average  value  of  the  terms  on  the  right  hand  side  of  (4-29) 
is  removed  at  each  time  step.   This  sets  the  constant  to 
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zero  and  does  not  allow  the  solution  to  grow  by  a  constant 
each  iteration. 

The  current  technology  of  large  mainframe  computers 
incorporates  vector  processing.   Utilization  of  this 
ability  can  be  achieved  when  the  system  of  equations  can  be 
vectorized.  The  GFEM  is  derived  in  vector  formulation  and 
is  easily  vectori  zed . 

D.    INITIAL  CONDITIONS 

1 .    Simple  Atmospheric  Waves 

The  initial  conditions  must  allow  a  relative 
amount  of  control  for  the  desired  input  parameters  as  well 
as  satisfy  the  boundcry  conditions.   The  forecast  model 
history-carrying  variables  are  (J),  u  and  v.   The  analytic 
expression  for  the  s  t  reanif  uncti  on,  i];,  is 


^    =  A  sin^  ^  sin  ^f^  -  U(y-y  .  .)  +  ^ 
W      L       -^  -^  m  1  d '   f 

where  A  =  amplitude  of  the  perturbation, 

W  =  width  of  the  channel, 

L  =  length  of  the  channel, 

n  =  wave  number, 

0  =  mean  flow  speed,  and 

^mid  "  middle  point  of  the  channel 


(4-30) 


40 


The  first  term  in  the  expression  is  the  perturbation. 
The  second  term  represents  the  north/south  slope  necessary 
to  support  a  mean  flow  of  0.  The  third  term  is  the  mean 
depth  term.   The  geopotential  height,  <^y    is  related  geostro- 
phically  to  the  streamfuncti  on ,  ^p  y    by 


(j)  =  fo'^ 


(4-31) 


By   use  of   a   trigonometric   identity,    the    streamfuncti on    is 
written    as 


A,,  27Ty,     .    2-nnx        .. ,  \    ^     * 

\b    =   ttCI    -    COS — ^)sin— i -    U(y-y    .  .)    +  :=- 

^2^  w    '  L  ^■^    •'mid'        f. 


(4-32) 


The   vorticity    is    given    by 


where 


C=    Za^^^A    sina2X    cos    2a]^y   -    a2^isina2X 
+   02       ^sina2X   cos   2a]^y 


0.2 


_   TT 


w 

27Tn 


(4-33) 


The  divergence  is  determined  through  a  linearized  form  of 

the  quasi -geostrophic  divergence  equation  from  Chapter  3-2 

of  Haltiner  and  Williams  (1980),  in  the  form 

f2      f 

/   °  =  /  U^^^  (4-34) 


V^D 
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The    right-hand    side   of   this    equation    becomes 

f  f  ^ 

cos   ot2X    [(-OflZai^Aae    +^Ua2^f)cos    Za^y    -    -yUae^l^  (4-35) 

Assuming  a  solution  form  for  divergence  on  the  left-hand 
side 

D  =  cos  agxCC^cos  2a]^y  +  C2) 

the  constants  C^^  and  C2  become 

f  .         f«     A  ^  ^ 

Ci  =  -  (-^U2a^2Ae,2  +  ^Ua23f )  /  ( 4ai2  +  ^^2^  +  -^  ) 


and 


f.   3  ,     2    f^" 
/^^2   f/(^2   ^^^ 


Using  v^x  =  D  and  assuming  a  solution  form  for  x  3s 


=  cosa2x( C3C0s2aiy  +  C4) 


(4-36) 


then  C3   and  C4   become 

C3  =  -  Ci/(4ai2  +  ^2^) 


and 


C4  =  -  02/02 
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Observe   that   i|^    equals    a    constant    and  -r^  equals    zero    along 
the   north/south   walls.      With   the   expressions    for    stream- 
function,     (4-32),  anc    velocity    potential,     (4-36),  appropri  ate 
formulae   for   j    and    v    can    be    derived. 

u   =    U    -    sin    a2x(Aai    sinZa^y   +   03(03    cos    Za^y   +    C4))     (4-37) 


and 


V  =  cos  a2x(a22(l-cos  Za^y)  -  2C3a^sin  2aiy)        (4-38) 

and  (j)  is  given  by  eiquations  (4-31)  and  (4-32). 
2.    Source  Term 

The  addition  of  a  source  term  into  the  continuity 
equation  will  test  the  resol vabi 1 i ty  of  the  GFEM  model  for  a 
source.   The  source  is  constructed  to  represent  a  small- 
scale  meteorol  cgi  ca';  phenomenon  to  provide  a  measure  of  the 
response  of  the  GFEfl  model  to  forcing  near  the  smallest  grid 
scale.   The  source  initial  conditions  must  be  able  to 
describe  a  small-scale  feature  embedded  in  a  large  synoptic 
situation.   The  initial  conditions  are  derived  by  choosing  a 
desired  solution  and  then  back  substituting  to  derive  the 
form  of  the  source  required.  This  is  a  unique  approach  when 
adding  a  source  term.   First,  a  source,  S,  is  added  to  the 
continuity  equation  in  the  form  of 


3t    8x    3y   ^^3x    3y^ 


<t> 


(4-39) 
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Following  the  scale  analysis  of  Chapter  3-2  in  Haltiner  and 
Williams  (1980), a  new  quasi -geostrophi c  potential  vorticity 
equation  is  formed  with  the  source  term  included, 


^Jt'  %^^%^)(?  -%  -  fs  = 


0 


Solving  for  S  leaves 


s  -  (J-      11  _L  +  li  J_wii      I/lliiL  4. 

^at    ay  3x    3x  3y' ^  $  "  f 4  2 

A 


LA)) 


(4-40) 


The  technique  for  testing  the  source  term  is  to  pick  a 
streamf uncti on  and  then  solve  for  the  source  term  which  will 
satisfy  that  streamf uncti on.   The  particular  streamf uncti on 
evaluated  will  be  a  steady  state  solution.   The  source  in 
thi  s  case  i  s 


2     2  2     2 

-    1  (d±     3  /  3  i>    ,    l_i\    3Jl  j_^l_i  +  l_jLn 


%''y    '^\x^    ■  3y2'    '^    ''\.^        3y2 


(4-41) 


The  particular  form  of  the  streamf uncti on  is 


*  =  -  u(y-y„,id'  -  '^=1"   ij^ 


1  .-     2,7TX, 


(4-42) 


where  R  is  a  constant.   The  use  of  the  exponential  allows  a 
variation  in  the  scale  whereby  a  small  scale  feature  can  be 
embedded  in  a  larger  flow  pattern.   Substitution  of  this 
particular  ^ii    into  the  expression  for  the  source  yields 
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1    /,",     j^    Ait     .     ZttV 


^   c  i  «  2      TT  X 


2        2  2  2 

f^^^^       c-Jn      ^y  /O 2ttX  le4«      2TrXx       « /.   /  TT  X  2  TTV  \ 

^rT^    Sin   -ijj^(2cos-^  -   ^sin  -i7-)-2A{^^)    cos-^) 


+   e 


-    R    '^"   —,A,.    2        2 


/A/TTx       ^,.  TTV,47r      .,-„   2ttX      ^      2  IT  .      4  TT  X  ,   w 

(rCl)    sin   -^C—  sin.-^-  +  _  sin-^})y 


1     C-;  „    TTX 
Sin 


-    f— (r-  l    ^^"    ^  W  ^^  sin -j— ) 


(e 


"I     ^-;«       TTX 
-R    '^" 


'^^•"^C^r^    (J)(2cos2lx. 


R    ^^"      ~L-^ 


+    4A(J)^)) 


■(4-43) 


E.    NONLINEAR  EFFECTS 

Cullen  (1974a)  discussed  the  requirement  to  sin.ulate 
the  cascade  of  energy  to  unresolved  wavelengths.   To  avoid 
accumulation  of  energy  in  the  shortest  resolvable  wave- 
lengths, he  developed  a  two-step  method  for  computing  the 
nonlinear  terms.   Implementation  of  the  two-step  method  can 
be  accomplisned  efficiently.   In  the  formal  development  of 
the  model  equations,  there  are    several  nonlinear  terms.   The 
two-step  method  involves  projecting  both  variables  in  a 
nonlinear  term  into  Galerkin  space  and  then  performing  the 
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multiplication.      This    procedure    insures    the   best    approxima- 
tion  to  that   nonlinear   term    in    a    weighted-mean    sense. 
Additionally,    the   computational    effort    for   the   two-step 
method   is   much   less  than   the   one-step  method   which    involves 
computation    of    all    the    possible    interaction    coefficients. 

F.    STABILITY  CRITERION 

The  stability  criterion  for  the  forecast  equations  is 


A  t  < 


1 


AX 


+  f 


This  is  more  restrictive  than  a  normal  Courant-Fredri ch-Lewy 
'CFL)  condition.   It  is  made  up  of  the  normal  CFL  criterion 
based  on  advection  plus  another  effeci'  due  to  rotation, 
(inertial  motion).   It  provides  an  upper  boundary  to  the 
maximum  allowable  time  step  when  using  a  semi- implicit 
scheme.   In  this  research,  time  steps  in  excess  of  an  hour 
are  easily  used.   Care  must  be  taken  not  to  exceed  the 
stability  criterion  when  larger  mean  flows  are  experienced 
or  when  moving  to  higher  latitudes  (larger  f ) . 

A  more  detailed  discussion  of  the  derivation  of  the 
stability  criterion  can  be  found  in  Appendix  A. 
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V.   EXPERIMENTS  AND  RESULTS 


Experiments  will  be  performed  involving  the  features  in 
Chapter  II  to  test  the  hypothesis.   Four  separate  experi- 
ments will  be  performed.   The  first  experiment  will 
establish  the  overall  performance  characteristics  of  three 
specified  models.   The  next  three  experiments  will  address 
each  specific  feature  supporting  the  hypothesis. 

Standards  of  comparison  must  be  established  for  each 

experiment.   The  basic  comparison  will  be  achieved  through 

harmonic  analysis  of  the  initial  and  forecast  fields.   The 

harmonic  analysis  is  a  double  Fourier  decomposition.   The 

harmonic  analysis  requires  that  the  data  be  on  a  uniform 

grid.   Because  of  the  triangular  and  variable  grids,  an 

interpolation  is  necessary.   A  fifth  degree  polynomial 

surface  fitting  routine  from  the  International  Mathematics 

and  Statistics  Library  (IMSL)  was  employed  to  interpolate 

data  on  the  non-uniform  grids  to  a  uniform  spacing.   First, 

a  harmonic  analysis  is  performed  along  columns  of  constant  x 

to  obtain  the  y  structure  of  the  initial  conditions.   Then  a 

harmonic  analysis  of  the  amplitudes  of  each  y  wave  number 

along  rows  gives  a  true  double  harmonic  analysis.   This 

double  transform  approach  is  necessary  due  to  the  initial 

2 
conditions.   The  use  of  sin  -g=^  in  the  y  direction  is  the 


2tt 


same    as    using    (1.0    -    cos    -^).      The    boundary    conditions    are 
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satisfied  by  a  sine  wave  in  the  y  direction.   Therefore,  an 
infinite  sine  series  is  needed  to  represent  the  cosine 
f uncti  on . 

The  theoretical  movement  of  each  wave  is  obtained  from 
an  analysis  of  the  shallow -water  equations.   The  analysis  in 
Section  2-6  of  Haltiner  and  Williams  (1980)  predicts  the 
wave  speed  as 


C  = 


Cf/H)|f 

v-  +  f^/gH 


(5-1) 


where 

C  equal s  wave  speed, 

H  equals  height,  and 

u  X  direction  wave  number. 

By   generalizing   the    streamf uncti on   to   a   two-dimensional    wave 
form 


li;    =    A    sin   — Jr^  cos    Ti(x-ct) 
w 


and  after  further  manipulation,  the  wave  speed  becomes 


C  =  Uei/Cl  +  ^o^/Hy/  +  Uy^))) 


(5-2) 


where 


$  equals  average  geopotential  height, 
yj^  equals  x  direction  wave  number,  and 
Uy    equals  y  direction  wave  number. 
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A.    EXPERIMENT  1 

Establishing  wave  propagation  characteristics  on  a 
uniform  grid  will  illuminate  the  properties  of  each  model. 
Three  models  will  be  compared: 

1.  A  triangular  subdivision  GFEM  model; 

2.  A  rectangular  subdivision  GFEM  model;  and 

3.  An  equivalent  finite  difference  model. 

The  GFEM  models  will  use  the  differentiated  form  of  the 
shallow-water  equations  with  a  semi -impl i ci t  scheme.   The 
finite  difference  model  uses  the  undifferentiated  form  on  a 
staggered  grid  with  a  centered  (leapfrog)  time  scheme.   The 
triangular  subdivision  uses  equilateral  triangles.   The 
relationship  between  base  and  height  for  an  equilateral 
tri  angl e  i  s 

BASE  =  2  X  HEIGHT//3 

The  X  distance  therefore  has  the  same  relationship  to  the  y 
d  i  stance 

AX  =  2  Ay//3 

To  perform  comparisons  of  similar  models,  the  same  x,  y 
relationship  is  retained  for  the  rectangular  and  finite 
di  f f erence  model s. 

The  basic  difference  between  the  rectangular  and 
triangular  models  will  be  in  the  approximating  polynomials. 
The  rectangular  polynomials  are   bilinear  while  the 
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triangular  polynomials  are  linear.   The  higher  order 
polynomials  should  provide  an  increase  in  accuracy.   During 
the  integration  process,  many  integrals  require  evaluation. 
Numerical  quadrature  could  be  utilized,  however,  a  more 
efficient  method  is  available  through  the  use  of  natural 
coordinates.  Quadrature  with  no  error  can  be  accomplished 
by  formula  with  this  method.   A  detailed  description  of  the 
natural  coordinate  method  is  delineated  in  Appendix  B  for 
both  triangles  and  rectangles. 

Diagrams  for  the  grids  for  the  triangular,  rectangular 
and  finite  difference  models  are  shown  in  Figures  2,  3  and  4 
respectively.   The  domain  is  4896.0  km  in  the  y  direction 
and  5653.0  km  in  the  x  direction.   Each  model  has  12  incre- 
ments in  the  x  and  y  directions.   This  gives  156  degrees  of 
freedom.  All  tests  will  be  for  a  48  hour  time  integration, 
a  mean  flow  of  10  m/s  and  a  mean  depth  of  1000  meters.  A 
perturbation  of  1  m/s  will  be  added  to  the  geopotential 
field,  which  includes  the  mean  height  plus  the  north/south 
slope  required  for  the  10  m/s  mean  flow.   This  small  pertur- 
bation will  focus  on  the  linear  aspects  of  the  model's 
capabilities.   Experiments  described  below  will  evaluate 
larger  perturbations,  and  hence  the  nonlinear  interactions. 
Models  are  evaluated  for  wave  numbers  1,  2,  3  and  4. 
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Fig.  3.    Rectangular  uniform  subdivision  for  a  channel  in 
Cartesian  coordinates. 
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Fig.  4.    Finite  difference  uniform  grid  for  a  channel  in 
Cartesian  coordinates. 
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B.    RESULTS  1 

The  propagation  characteristics  of  the  GFEM  models  and 
the  finite  difference  model  on  a  uniform  grid  will  be 
evaluated.   Graphical  displays  will  highlight  synoptic 
characteristics.   Harmonic  analysis  will  show  phase 
propagation  cind  wave  amplitude  distributions. 

The  initial  and  48-h  forecasts  for  the  conditions 
stated  in  experiment  1  with  wave  number  one  on  a  triangular 
grid  are    shov»n  in  Figures  5  and  6.   The  southeast/northwest 
orientation  in  the  forecast  fields  is  due  to  the  sine  wave- 
form in  the  y  direction.   This  waveform  with  model  boundary 
conditions  requires  an  infinite  number  of  y  modes. 
Fortunately,  the  amplitudes  of  the  higher  wave  numbers  are 
small.   Consequently,  there  are    only  a  few  important  modes. 
The  multiple  modes  cause  a  skewness  in  the  solutions  for  all 
three  models  since  the  y-modes  have  different  phase  speeds. 
Additionally,  an  analytic  expression  for  the  solution 
requires  an  infinite  sum  of  the  modes.   For  this  reason,  no 
true  solution  has  been  computed  to  compare  with  the  GFEM  and 
FDM  model  s. 

The  divergence  equation  is  the  most  sensitive  equation 
of  the  three  forecast  equations.   A  mean  depth  of  1000 
meters  was  chosen  because  the  divergence  is  large.   There- 
fore, a  critical  evaluation  of  the  sensitivity  of  the  GFEM 
model  to  this  important  meteorological  parameter  can  be 
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Fig.  5.    Initial  conditions  for  the  GFEM  model  with  a 
triangular  subdivision  and  wave  number  one. 
Contour  intervals  are  600  m^/s^  for  geppotential 
height,  .2  m/s  for  u  and  v,  .6x10"°  s"^  for 
vorticity  and  .2xlO-7s-l  for  divergence.   Nodal 
points  are  denoted  by  an  x. 
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Fig-    6.         As    in    Fig.    5    but    a   48-h    forecast 
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performed.   A  divergence  field  which  has  propagated 
uniformly  with  the  vorticity  field  is  shown  in  Figure  6. 
Harmonic  analysis  (not  shown)  indicates  divergence  ampli- 
tudes which  initially  adjust  and  then  oscillate  slightly 
about  a  mean  value  of  .1270x10"'  s''^.   The  phases  of  all 
forecast  fields  are  uniform.   Synopti cal ly ,  the  forecast 
fields  4),  u,  and  v  indicate  a  smoothly  resolved  atmospheric 
wave. 

A  rectangular  subdivision  48-h  forecast  for  wave  number 
one  is  shown  in  Figure  7.   Comparisons  between  Figures  6  and 
7  show  little  if  any  difference  which  indicates  that  the 
triangular  and  rectangular  subdivisions  give  comparable 
results.   Harmonic  analysis  of  the  divergence  fields 
indicates  a  small  oscillation  about  a  mean  similar  to  that 
observed  for  the  triangular  subdivision.   The  48-h  forecast 
from  the  finite  difference  model  for  wave  number  one  is 
shown  in  Figure  8.   The  equivalent  finite  difference  model 
also  has  12  increments  in  each  direction.   Comparisons 
between  Figures  6,  7  and  8  show  that  the  finite  difference 
model  vorticity  field  lags  the  GFEM  models.   Additionally, 
the  divergence  field  has  lost  its  areal  extent  and  developed 
a  strong  gradient  of  divergence. 

In  each  of  the  wave  number  one  tests  there  were  12  grid 
points  representing  the  wave.   Good  propagation  should  be 
expected  from  all  three  of  the  models  for  a  12-increment 
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Fig.  7.    As  in  Fig.  6  but  a  rectangular  subdivision 
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Fig.    8.         As    in    Fig.    6    but    a    FDM   model 
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wave.   As  the  wave  number  is  increased  with  the  same  degrees 
of  freedom,  the  wave  resolution  will  be  decreased,  and 
clearer  comparisons  can  be  made. 

Figures  9  and  10  are  the  initial  and  48-h  GFEM  model 
forecasts  for  wave  number  2,  triangular  subdivision. 
Forecasts  for  the  rectangular  subdivision  and  the  FDM  model 
are  shown  in  Figures  11  and  1?..      The  triangular  and 
rectangular  subdivision  forecasts  are  similar.   However,  the 
FDM  model  forecast  displays  additional  divergence  and 
vorticity  centers  near  the  boundaries.   Figures  13  through 
16  are   a  series  similar  to  Figures  9  through  12  except  for 
wave  number  three.   These  waves  are  represented  by  4  grid 
points.   There  are  indications  in  Figure  14  that  energy  is 
appearing  in  wave  numbers  other  than  wave  number  3. 
Comparison  of  Figures  14  and  15  show  that  the  triangular 
forecast  contains  more  noise  than  the  rectangular  forecast. 
Figure  16  shows  a  divergence  field  from  the  FDM  model  which 
is  very  different  from  Figures  14  and  15.   There  is  a 
north/south  asymmetry  in  the  u  and  divergence  fields.   The  u 
cells  in  the  northern  regions  appear  to  be  expanding  while 
the  southern  cells  are  decreasing. 

Table  1  is  a  composite  of  the  harmonic  analysis  for 
the  triangular,  rectangular  and  finite  difference  wave 
number  3  case  (v  component  amplitudes).   The  initial  and 
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Fig.    9.        As    in    Fig.    5    except    for    wave    number    two.      Contour 
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Fig.    10.      As    in    Fig.    9    but    a   48-h    forecast. 
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Fig.  11.   As  in  Fig.  10  but  a  rectangular  subdivision 
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Fig.    12.      As    in    Fig.    10    but    a    FDM   model 
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Fig.  13.   As  in  Fig.  5  except  for  wave  number  three. 

Contour  intervals  are  600  m^/s^  for  geopotential 
height,  .4  m/s  for  u,  2.0  m/s  for  v,  .6x10-5  s-1 
for  vorticity  and  .2x10"^  s'^  for  divergence. 
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Fig.    14.      As    in    Fig.    13    but    a   48-h    forecast 
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Fig.    15.      As    in    Fig.    14    but    a    rectangular    subdivision 
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Fig.    16.      As    in    Fig.    14    but    a    FDM   model 
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48-h  V  component  amplitudes  for  wave  number  3  in  the  x 
direction  and  the  first  six  possible  y  modes  are  tabulated 
bel ow. 

TABLE  1 
Harmonic  analysis  (v  component  m/s)  -  wave  number  3 


Tri  ang 

les 

Wave 

Ini  ti  al 

48-H 

1 

2.12 

2.14 

2 

.01 

.04 

3 

.41 

.43 

4 

.00 

.00 

5 

.05 

.06 

6 

.00 

.00 

Rectangl es 
Initial    48-H 


2.20 

2.22 

.01 

.02 

.44 

.44 

.00 

.00 

.06 

.06 

.00 

.00 

Finite  Difference 
Initial    43 -H 


2.20 

2.20 

.01 

.09 

.44 

.42 

.00 

.00 

.06 

.05 

.00 

.00 

The  minor  differences  in  the  initial  state  amplitudes 
are  due  to  the  differences  between  the  triangular  and 
rectangular  approximation  for  the  initial  conditions. 
However,  the  48-h  amplitudes  show  that  the  rectangular 
version  has  slightly  better  preserved  the  initial  state 
amplitudes  than  either  the  triangular  or  finite  difference 
models.   The  finite  difference  model  has  transferred  more 
amplitude  to  the  other  modes,  especially  in  the  second  y 
mode. 

Figures  17  and  18  are  the  initial  and  48-h  GFEM  model 
triangular  subdivision  forecasts  for  wave  number  4.   Wave 
number  4  is  a  severe  test  for  this  model.   The  energy 
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Fig.    17.      As    in    Fig.    13    except    for   wave    number    four 
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Fig.   18.      As    in    Fig.   17    but    a   48-h    forecast. 
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transfer  to  other  modes  is  readily  apparent.   Additionally, 
both  the  u  component  of  the  wind  and  the  divergence  show 
considerable  noise  at  the  boundaries.   Figure  19  is  the  48-h 
forecast  for  the  rectangular  model.   No  energy  transfer  nor 
noise  at  the  boundaries  is  apparent  in  Figure  19.   Figure  20 
is  the  48-h  forecasts  for  the  finite  difference  model.   The 
asymmetry  in  the  wave  number  three  finite  difference  case 
is  magnified  for  wave  number  four.   The  divergence  field  is 
heavily  asymmetric.   The  u  cells  have  developed  a  dipole 
from  each  original  center.   Comparisons  of  Figures  18,  19 
and  20  show  the  superiority  of  the  rectangular  forecast. 
While  the  triangular  version  is  noisy,  it  still  has 
maintained  the  main  features  better  than  the  finite 
di  f f erence  model  . 

Table  2  is  similar  to  Table  1  except  for  the  wave  number 
4  case. 

TABLE  2 
Harmonic  analysis  (v  component  m/s)  -  wave  number  4 

Triangles  Rectangles  Finite  Difference 
Ini  ti  al 
2.94 
.01 
.58 
.00 
.08 
.00 


Wave 

Ini  ti  al 

48-H 

1 

2.68 

2.60 

2 

.00 

.03 

3 

.51 

.59 

4 

.00 

.00 

5 

.06 

.08 

6 

.00 

.00 

48-H 

Ini  ti  al 

48-H 

2.96 

2.94 

2.96 

.10 

.01 

.05 

.56 

.58 

.60 

.04 

.00 

.05 

.08 

.08 

.08 

.00 

.00 

.01 
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Fig.  19.   As  in  Fig.  18  but  a  rectangular  subdivision 


72 


GEOPOTENTIRL  HEIGHT 


STREAMLINES 


X-flXIS 

V 

TflU  -  i9 


.^„f- — — /"\  '      '/ 

!«"\x     x/x\x     x/x\x     xj   X  \x     xf 


2.3  3.0  4.0 

X-flXIS    IfCTERS) 


0.0  IX  3.0  3.0  4.0  5.0 

X-flXlS  (METERS)  hIQ* 


VGRTICITY 


DIVERGENCE 


3.0  3.0  4.0 

X-flXIS   I METERS) 


0.0  1^  3.0  3.0  4.0  5.0 

X-HXIS   t METERS)  hIQ* 


Fig.    20.      As    in    Fig.    18    but    a    FDM   model 
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The  harmonic  analysis  shows  that  the  triangular  version 
transfers  more  energy  in  the  fifth  y  mode  then  does  either 
other  model.   However,  the  rectangular  model  transfers  more 
energy  in  the  second  y  mode  than  does  either  other  model. 
Table  3  compares  the  divergence  for  the  wave  number  4  case. 


TABLE  3 


Harmonic    analysis    (divergence    s"-*-)    -    wave   number   4 
(All    numbers    are    scaled    10    to   the   minus    10.) 


Wave 
1 
2 
3 
4 
5 
6 


Tri  angl es 

Initial         48-H 


2682.0 

.0 

521.3 

.0 

62.5 

.0 


391.1 
20.5 
60.5 
10.5 


Rectangles  Finite   Difference 

Initial         48-H         Initial         48-H 


1741.0   2929.0 
140.8       .0 


589.8 

.0 

84.4 

.0 


2179.0   2801.0 
207.5    447.4 


440.7 

81,2 

56.8 

8.4 


545.5 

.0 

68.7 

.0 


4462.0 
5934.0 
2314.0 

923.5 
1382.0 

882.0 


The    finite   difference   model    has   greatly    amplified   the 
initial    divergence.      A   choice   between    the   triangular   or 
rectangular    grids    version    based    only    on    these    results    would 
be    purely    subjective.      However,    the    response   of   both    grids 
at    this    high    wave    number   in   terms   of  divergence   is    highly 
superior   to   the    FDM. 
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Figure  21  is  the  phase  propagation  diagram  for  the 
three  models  obtained  from  the  double  harmonic  analysis.   In 
all  cases,  the  percentage  propagation  is  for  the  first  y 
mode  and  the  appropriate  x  mode.  The  phase  propagation 
diagram  indicates  that  the  triangular  and  rectangular  models 
are  con-parable  for  the  low  wave  numbers  and  diverge  slightly 
at  higher  wave  numbers.   In  all  cases  the  GFEM  models  have 
better  wave  propagation  than  the  FDM  model. 

The  analysis  performed  on  the  uniform  grids  has  shown 
that  each  GFEM  model  performs  in  a  highly  satisfactory 
fashion.   There  are  minor  differences  in  the  phase 
propagation,  while  the  rectangular  version  has  an  advantage 
in  the  control  of  energy  transfers  which  manifest  themselves 
in  the  divergence.   The  finite  difference  model  performs  as 
well  as  the  GFEM  models  for  the  long  waves.   However,  when 
forecasting  the  shorter  waves,  the  finite  difference  model 
does  not  perform  as  well  as  either  GFEM  model. 

C.    EXPERIMENT  2 

This  experiment  will  investigate  two  options  for  a 
variable  grid.   Atmospheric  wave  propagation  on  variable 
triangular  and  rectangular  grids  will  be  the  test  vehicle. 
Since  this  FDM  model  has  no  capability  on  a  variable  grid, 
it  will  not  be  evaluated  during  this  experiment.   The 
comparisons  will  be  between  a  rectangular  and  triangular 
GFEM  model  . 
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Fig.  21.  Phase  propagation  diagram  for  the  triangular, 
rectangular  and  finite  difference  models  from 
Experiment  1. 
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There  are  many  criteria  for  determination  of  a  suitable 
variable  grid.   Both  rectangles  and  triangles  allow  a 
successful  implementation  of  a  variable  grid.  One 
criterion  is  the  accuracy  of  resolving  the  atmospheric  waves 
that  move  through  the  variable  grid.   A  second  criterion  is 
the  increase  in  resolution.   Another  consideration  is  the 
ease  with  which  the  grid  can  be  refined. 

Cullen  (1979)  stressed  the  importance  of  a  smooth 
variation  in  element  size  when  increasing  the  resolution 
from  coarse  to  fine.   Older  (1981)  developed  a  technique  for 
transforming  a  uniform  triangular  grid  into  a  variable  grid 
when  working  with  a  GFEM  model.   He  demonstrated  that  a 
smoothly  varyirg  grid  allowed  atmospheric  waves  to  propagate 
through  the  transition  zone  with  much  less  noise  generation 
than  an  abruptly  varying  grid.   Some  finite  difference 
models  with  nested  grids,  for  example  Harrison  (1973),  have 
abruptly  varying  resolution.   A  generalization  of  Older's 
techniques  allows  a  similar  transformation  for  a  rectangular 
grid. 

The  test  cases  for  experiment  2  will  be  for  a  simple 
atmospheric  wave.  A  mean  depth  of  5  km,  a  mean  flow  of 
10  m/s,  wave  number  one  and  a  perturbation  of  1.0  m/s  are 
empl oyed . 
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D.    RESULTS  2 

The  ability  to  refine  a  grid  should  be  easily  achieved, 
The  method  of  transforming  from  a  uniform  grid  to  smoothly 
varying  should  allow  one  to  choose  not  only  what  degree  of 
refinement  but  also  where  the  refinement  occurs.   Older 
(1981)  developed  a  transformation  method  based  on 

X  =  X  +  A  cos  Kx 


where 

X  =  the  transformed  grid 

X  =  the  original  grid 

A  =  a  constant 

k   =    Ztt/L. 
Woodward    (1981)    modified    the   transformation    to    include    a 


SI  n 


kx 


for  the  longitudinal  stretching.   A  trignometric  identity 
is  empl oyed  yielding 

X  =  X  +  A(l  -  cos  kx) 
This  research  added  the  capability  to  control  the  location 
of  the  refinement  through 

X  =  X  +  A(l  -  cos(kx  +  5 )) 

sX 
The   map   factor  r—   is   defined    by 

a  X 


3X 
3X 


=  1  +  kA  sinCkx  +  6) 
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The  maximum  and  minimum  values  are 


max  mm 

The  ratio  R  of  maximum  map  factor  to  minimum  map  factor  is 


R  = 


T+kA 
1-kA 


Sol vi  ng    for   A  yi  el ds 


A   = 


R-1 
k(R+l) 


For    purposes    of   this    experiment   R   was   chosen   to   vary   from 
1.0  to  4.0.      A  value  of  R  ^-  4   implies  that  the  minimum   X   is 
one    fourth    of   the   maximun;i   X.      The   Y   transformation    is 
performed  in  a  similar  fashion  except  that  Y  =  y  +  B   sin  Ly 
is  empl oyed   where 

Y  =  the  transformed  grid, 

y   =    the   original    grid, 

B   =    a    constant,    and 

L   =    2Tr/W 
Placement    of   the    high    resolution    is   accomplished   by 
determining  the   value   of  6    required   to   place   the   minimum   map 
factor    as   desired.      For    instance,    if   the    refined    grid    is 
desired    in    the   middle   of   the   channel    (L/2)    then   <5    would    be 
equal    to   tt/2. 
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Figure  22  is  the  48-h  forecast  with  the  triangular 
subdivision  for  a  uniform  grid.   This  forecast  will  act  as 
the  control  for  the  triangular  cases.   Figure  23  is  the 
triangular  subdivision,  48-h  forecast  with  R  =  2.0.   Figure 
24  is  the  triangular  subdivision,  48-h  forecast  with 
R  =  3.0.   Figure  25  is  the  triangular  subdivision,  48-h 
forecast  with  R  =  4.0.   All  of  the  48-h  u  and  v  fields 
appear  identical.   This  indicates  that  there  is  no  degra- 
dation in  forecast  skill  when  using  the  same  number  of 
degrees  of  freedom  and  locating  many  of  those  unknowns  into 
a  refined  area.   There  is  no  noise  apparent   in  the 
transition  regime.   However,  there  are  definite  differences 
in  the  divergence  field,  although  the  magnitude  of  the 
maximum/  minimum  value  of  divergence  is  relatively  small 
(1.0x10'^  s'b. 

Table  4  shows  the  harmonic  analysis  amplitudes  of  the 
geopotential  fields  for  the  different  triangular  cases.   The 
relationship  between  initial  and  final  amplitudes  remains 
the  same  regardless  of  the  degree  of  variability.   The  phase 
propagation  of  the  v  field  in  the  control  case  is  97.7%, 
whereas  it  is  96.7%,  95.4%  and  94.1%  for  R  values  of  2,  3 
and  4.   The  variation  from  the  control  of  the  phase  propaga- 
tion with  grid  refinement  is  under  four  percent. 
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Fig.  22.   As  in  Fig.  6  for  R  =  1.0.   Contour  intervals  are 
600  m^/s^  for  geopotential  height,  .2  m/s  for  u 
and  V,  .6x10"°  s"-^  for  vorticity  and  .6xl0"^s"^ 
for  di  vergence . 
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Fig.    23.      As    in    Fig.    22    for    R    =    2.0 
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Fig.    24.      As    in    Fig.    22    for    R    =    3.0. 
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Fig.    25.      As    in    Fig.    22    for    R    =    4.0. 
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TABLE    4 
Harmonic   Analysis    (Geopotenti al    m^/s^) 
Tri  angl es 
R    =    1.0  R    =    2.0  R    =    3.0  R    =    4.0 

Wave    Initial    48-H      Initial    48-H      Initial    48-H      Initial    48-H 
1         67.98      67.16      67.91      64.65      67.91      64.95      67.91      63.70 

.69  .00 


.00 


55 


.00 


.99 


.00 


.00  .21  .00  .37  .00 

1.90        1.70        1.67        1.30        1.33 

.00  .16  .00  .23  .00 


.57 


13.55      12.90      13.45      12.55      13.41      13.24      13.29      11.16 


.54  .00  .28 

1.62        1.32        1.38 

.36  .00  .21 


Figure  26  is  the  control  case  for  the  rectangular 
version,  48-h  forecast.   Figure  27  is  thv^  rectangular 
subdivision,  48-h  forecast  with  R  =  2.0.   Figure  28  is  the 
rectangular  subdivision,  48-h  forecast  with  R  =  3.0. 
Figure  29  is  the  rectangular  subdivision,  48-h  forecast 

with  R  =  4.0.   Close  inspection  reveals  no  noise  in  either 
the  <t>i    u  or  V  fields.   No  degradation  in  forecast  skill  has 
been  observed  due  to  increased  resolution  with  a  rectangular 
subdivision.   All  <|),  u  and  v  fields  are  very  similar  in 
structure.   The  main  difference  in  the  forecasts  lies  in  the 
divergence  field  where  the  magnitudes  are  small  (  10"°s"M 
because  of  the  specified  mean  depth. 
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Fig.    27.      As    in    Fig.    26    for    R    =    2.0 
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Fig.    28.      As    in    Fig.    26    for    R    =    3.0 
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Table  5  is  similar  to  Table  4  except  for  the 
rectangular   cases.      Phase    propagation    for   the   control    case 
V   field   is   98.0%,   whereas   it   is  98.3%,   98.1%,   and  97.9%   for 
R   values   of   2,    3,    and   4.      The   variation    from    the   control 
case   of   the    phase    propagation    with   grid   refinement    is    less 
than    half    a    percent. 

TABLE    5 
Harmonic   Analysis    (Geopotenti al    m^/s^) 
Rectangl es 
R    =    1.0  R    =    2.0  R    =    3.0  R    =    4.0 

Wave    Initial    48-H      Initial    48-H      Initial    48-H      Initial    48-H 
1         68.03      67.48      67.98      65.26      68.22      66.45      68.14      66.34 
.01         1.00  .03  .36 


.00 


.85 


.00        1.22 


13.59      12.95      13.39      13.26      13.49      13.81      13.36      12.39 


.01  .46  .05  .22 

1.92        1.81        1.57        1.64 

.02  .23  .07  .04 


.01  .48  .00  .43 

1.54        1.59        1.46        1.19 

.01  .30  .00  .22 


The  tests  thus  far  indicate  a  successful  grid  refine- 
ment capability  for  both  triangles  and  rectangles.   There 
are  neither  drastic  impairments  nor  improvements  in  either 
phase  propagation  or  amplitude  changes  for  either  technique. 
Neither  version  stands  out  above  the  other  although  the 
rectangular  version  has  less  phase  propagation  variation  and 
better  amplitude  conservation.   Both  versions  allow  compara- 
ble grid  refinement  through  the  use  of  the  procedure 
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previously   described.      Both    allow    an    easy    and    straightfor- 
ward  implementation   of   the    procedure. 

However,    there   is   a   more   fundamental    difference   between 
the   versions.      Operators,    such    as    a    simple   derivative,    have 
weighting   coefficients   associated   with    finite   element    appli- 
cations  as   well    as   with    finite   difference    scheme:;.      For 
instance,    a    simple-centered,    one-dimensional    finte   differ- 
ence  derivative   has    weights    of   1/2    and    -1/2.      The    finite 
element    weights    are    identical    to   the   finite   difftjrence 
weights   in   this   one-dimensional    case.      The    weights    can 
similarly   be   computed   for   two-dimensional    operatDrs    such   as 
a   Laplacian.      The   quadrature    formulas    used   to   determine   the 
weights    are    explained    in    more   depth    in    Appersdix  I).      The 
geometry   of   the   triangles    affects    the    weights.      Tor 
instance,    the   two-dimensional    Laplacian    applied    to    a   finite 
element   application   on    an   equilateral    triangular    subdivision 
gives    weights    of    1/6    to    all    surrounding    points,    and    -1    to 
the   center    point,    as    indicated    in    Figure   30.      However,    the 
weights  change  if  a  triangle  with  a  height  equal   to  one  half 
the    base    is    employed    as    shown    in    Figure   31.      The    fact    that 
some  weights  are  zero  is   a   strong  cause  for  alarm.      If  the 
triangles    are    flattened    further,    the    farthest    end    point 
weights    become    negative.      This    is    strictly    a    function    of   the 
geometry  of  the  triangle  and  is  a  clear  sign  that  use  of 
the    triangles    warrants   extreme   caution.      In    fact,    Williams 
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Fig.  30.   Woiighxs  associated  with  a  Laplacian  operator 
for  a  triangular  subdivision  using  equilateral 
tri  angl es. 


Fig.  31.  Weights  associated  with  a  Laplacian  operator  for  a 
triangular  subdivision  using  triangles  with  a  base 
equal  to  twice  the  height. 
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(unpublished  notes)  has  determined  that  the  geometry  of  the 
triangles  dictates  that  unless  a  triangle  with  at  least  t»vo 
equal  sides  (isosceles)  is  employed  the  boundary  conditions 
will  not  be  exactly  satisfied. 

E.    EXPERIMENT  3 

The  capability  of  the  GFEM  model  to  resolve  atmospheric 
phenomena  near  the  scale  of  the  smallest  grid  length  will  be 
demonstrated.   This  experiment  is  in  response  to  Task  II  of 
the  hypothesis. 

The  ultimate  motivation  for  any  "improvement"  in  a 
model  is  a  better  forecast.   During  the  past  decade, 
improvements  have  come,  to  a  certain  degree,  by  increased 
resolution.   While  it  is  intuitively  obvious  that  increased 
resolution  will  improve  a  forecast,  demonstrations  have  not 
been  made  to  show  how  the  model  resolves  features  near   the 
smallest  grid  length.   Staniforth  and  Mitchell  (1978) 
increased  the  resolution  with  a  variable  grid  but  resolved  a 
synoptic  scale  feature  in  a  37x37  uniform  fine  mesh  area. 

In  this  experiment,  a  source  term  simulating  a  mass 
source  is  added  to  the  continuity  equation.   An  expression 
for  the  source  term  in  terms  of  a  wave  form  is  derived  by 
theoretical  means.   The  wave  form  allows  a  very  small  scale 
feature  to  be  simulated.   A  detailed  description  of  this 
procedure  can  be  found  in  Chapter  IV-D-2.   By  assuming  the 
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wave  is  in  equilibrium  and  steady  state,  the  solution  and 
source  are  known  for  all  time  when  the  Rossby  number  is 
smal  1  . 

The  development  of  the  source  term  expression  followed 
traditional  quasi -geostrophi c  theory.   One  requirement  for 
quasi -geostrophi c  theory  is  that  the  Rossby  number  be  small. 
In  the  linear  case,  the  Rossby  number  will  be  approximately 
0.1.   However,  in  the  nonlinear  case  the  Rossby  number  is 
approximately  0.4  and  the  quasi -geostrophi c  theory 
assumption  is  violated. 

The  first,  test  will  be  on  a  uniform  grid  for  the 
rectangular  GFEM  aid  finite  difference  models  and  the  second 
test  will  be  on  a  variable  rectangular  grid  GFEM  model. 
High  resolution  versions  of  the  uniform  test  will  be  run  and 
act  as  the  control.   A  mean  depth  of  5  km  was  selected. 
Perturbations  of  2.5  m/s  and  25  m/s  upon  a  mean  flow  of 
10  m/s  will  be  evaluated  for  the  uniform  grid.   This  will 
allow  both  the  linear  and  nonlinear  aspects  of  the  different 
models  to  be  observed.   Here  linearity  is  implied  by  the 
smallness  of  the  perturbation  amplitude.   The  source  term 
solution  is  nonlinear  and  nonlinear  interactions  do  occur 
when  the  amplitude  is  small.   Each  test  will  be  integrated 
for  96  h  of  forecast  time.   A  perturbation  of  25  m/s  upon  a 
mean  flow  of  10  m/s  will  also  be  evaluated  for  the  variable 
grid. 
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The  main  thrust  of  this  experiment  is  to  determine  the 
grid-scale  sensitivity  for  the  GFEM  model  and  the  equivalent 
finite  difference  model.   Sensitivity  will  be  measured  by 
insertion  of  a  source  term  into  the  continuity  equation. 
However,  the  theoretical  development  presented  in  Chapter 
IV-D-2  is  for  the  rotational  component  of  the  wind.   Because 
the  initial  conditions  for  u  and  v  are  non-divergent,  the 
divergence  in  the  model  must  increase  from  a  zero  initial 
state.   This  increase  will  require  an  adjustment  process. 
This  serendipity  effect  will  be  exploited  during  the 
compari  sons. 

The  linear  cases  examined  represent  a  small-scale  short 
wave  embedded  in  a  long  wave  pattern.  Using  a  small  per- 
turbation for  the  first  tests  will  establish  a  level  of 
confidence  that  each  model  responds  to  the  source  term 
creditably.   The  nonlinear  cases  represent  an  active  vortex 
in  a  mean  flow.  The  maximum  u  component  is  35  m/s  which  is 
hurricane/typhoon  strength.   Physically,  the  source  term  in 
the  divergence  field  opposes  the  advection  of  the  vorticity 
field,  and  "anchors"  the  steady-state  vorticity  solution. 

In  either  the  linear  or  nonlinear  tests,  the  source 
term  is  nonlinear.   By  constraining  the  perturbation  to  be 
small,  the  source  term  plays  a  linear  role.   However,  when 
the  perturbation  is  large,  the  source  term  allows  full 
nonlinear  interactions. 
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F.    RESULTS  3 

1 .    Linear  Case  -  Uniform  Grids 

Figures  32  and  33  are  the  initial  and  96-h  fore- 
casts for  the  rectangular  GFEM  model.   There  are  156  degrees 
of  freedom  in  these  forecasts.   The  initial  v  field  has  a 
maximum  value  of  +2.27  m/s  acros:;  a  separation  of  4  incre- 
ments.  The  initial  u  field  maximum  and  minimum  are  12.5  m/s 
and  7.5  m/s.   The  source  tern  is  shown  in  the  initial  fields 
and  remains  constant  throughout  the  integration.   The  diver- 
gence value  increased  to  a  steady  state  solution  during  the 
first  12  hours  with  an  oscillation  that  died  out  after  hour 
24.   Figure  34  is  the  equivalent  finite  difference  96-h 
forecast.   Figure  34  demonstrates  that  the  finite  difference 
model  performs  well  for  these  initial  conditions.   Figure  35 
is  the  high-resolution  96-h  forecast  for  the  rectangular 
GFEM  model.   Close  inspection  of  Figures  33  and  35  shows 
little  difference.   The  low  resolution  forecast  with  156 
degrees  of  freedom  has  converged  to  the  high  resolution 
forecast  with  600  degrees  of  freedom.   Figure  36  is  a  graph 
of  V  component  amplitudes  at  hour  96  as  a  function  of  x  and 
y  wave  numbers  for  the  low  and  high  resolution  and  finite 
difference  models.   The  high   resolution  (S253115)  and  low 
resolution  (S153115)  forecasts  are  in  excellent  agreement. 
However,  the  FOM  (F153115)  forecast  departs  from  the  control 
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Fig.  32.   Initial  conditions  for  the  GFEM  model  with  a 

rectangular  subdivision  and  source  term  added  to 
the  continuity  equation.   Resolution  is  low  and 
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even    with    this    relatively    linear   case.      For    uniformity    of 
comparison,    similar   ranges   of  x,   y   wave   numbers   will    be 
shown    in    Figures   36    and   41. 

2.        Nonlinear   Case   -    Variable   Grids 

Figures   37    and   38    are   the    initial    and   96-h 
forecasts  for   the   rectangular  GFEM   model    with   low 
resolution.      The    initial    v    field    has    a   maximum    of   +22.7    iri/s 
across  a  separation  of  4  increments,  but  most  of  the  change 
occurs    over   2    increments.      The   initial    u   field    maximum/ 
minimum  is  35  m/s  and  -15  m/s.     The  source  term  as  shown 
remains   constant    during    the    integrations.      The    adjustment 
process    is    evident    when    viewirg   96-h    forecasts.      The    final    v 
field    has    a    maximum    of   44.3    m/s,    and    a   minimum   of   -38.0    m/s. 
The  final   u  field  has  a  maximum  of  40.2  m/s  and  a  minimum  of 
-29.8    m/s.       Figure    39    is    the    equivalent    finite    difference 
model    forecast   of   Figure  38    ami    the   marked   differences    in 
the    forecast    are    readily   .apparent.      The   maximum    final    v 
component    is   41.2    m/s    and    the   minimum    is    -47.3    m/s.      The 
maximum   final    u   component   is   50.6   m/s   and   the  minimum   value 
is    -38.6    m/s.      Figure   40    is    the   96-h    forecast    for    the 
rectangular   high-resolution   GFEM   model.      Comparisons   of 
Figures  38,  39  and  40  show  a  convergence  of  the  low- 
resolution   GFEM    version   towards   the   high-resolution    solution 
with   a   much   poorer   showing    for   the    finite   difference   model. 
Figure  41   is  a  graph  of  v  component  amplitudes  at  hour  96  as 
a  function  of  x  and  y  wave  numbers  for  the  low  and  high 
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104 


GEOPOTENTinL  HEIGHT 


STREAMLINES 


0.0         1.0  3.0  3.0  i.a  s.o 

X-flXIS  (rtlERS)  xlO" 


X-flXIS 


u 


inu  -  96 
rsisaus 


T  "' ^___js^  X   J(.   »■'  xV-»i;^_Y     X i-     •'-' 


0.0  1.0  3.0  3.0  4.C  S.O 

X-RXIS    (METERS)  -lO* 


0.0  1.0  3.0  3.0  4.0 

X-FIXI5  (MEIERS) 


s.o 

HlQ' 


VGRTICITY 

TRU  -  95 

Fsisaus 


DIVERGENCE 


Q.O  1.0  3.0  3.0  4.0  5.0 

X-flXIS   (METERS)  xltf 


3.0  3.0  4.0 

X-FIXIS  (METERS) 


Fig.    39.      As    in    Fig.    38    for    a    FDM   model 
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Fig.  41.   V  component  amplitude  as  a  function  of  x  and  y 

wave  number  at  hour  96  for  the  high  (S258115)  and 
low  (R158115)  resolution  GFEM  models  and  the 
finite  difference  (F158115)  model.   Perturba- 
tion =  25.0  m/s.   Contour  interval  is  0.5  m/s. 
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resolution  and  finite  difference  model.   As  in  the  linear 
case,  there  is  strong  agreement  between  the  low  and  high 
resolution  forecasts.   The  departure  of  the  FDM  from  the 
control  is  readily  apparent  in  this  highly  active  case.   The 
inability  of  the  FDM  to  properly  handle  the  nonlinear 
interactions  manifest  itself  as  large  spurious  amplitudes  at 
high  X  and  y  wave  numbers. 

The  relatively  poor  showing  of  the  finite  difference 
model  for  the  forced  cases  shown  so  far  can  be  improved,  to 
a  certain  degree,  by  increasing  the  resolution.   However, 
increased  resolution  requires  more  computational  effort. 
Figure  42  is  the  96-h  forecast  with  a  perturbation  of  25.0 
m/s  for  the  FDM  model.   The  initial  conditions  are  identical 
to  the  nonlinear  cases  shown  in  Figure  37.   The  difference 
in  the  forecasts  is  the  resolution  where  there  are  now  576 
degrees  of  freedom  (24X24).   The  96-h  forecast  is  certainly 
better  than  the  12X12  forecast.   However,  there  is  some  high 
frequency  noise  in  the  forecast.   Figure  43  is  identical  to 
Figure  42  except  that  there  are  1296  degrees  of  freedom 
(36X36).   The  high  frequency  noise  is  readily  apparent  and 
is  a  manifestation  of  nonlinear  aliasing.   For  this  forced 
case  the  finite  difference  model  is  unable  to  achieve  the 
same  forecast  as  the  low  resolution  GFEM  model.   This  result 
will  be  amplified  when  a  model  with  the  same  number  of 
degrees  of  freedom  as  the  low  resolution  model  is  employed 
but  with  a  variable  grid. 
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Fig.    42.      As    in    Fig.    39    with    576    degrees    of    freedom 
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This  comparison  of  the  finite  element  and  finite 
difference  methods  for  the  source  term  dramatically  high- 
lights the  potential  improvement  when  utilizing  a  Galerkin 
approach.   Both  graphical  displays  and  harmonic  analysis 
leave  little  doubt  as  to  the  superiority  of  the  GFEM  model 
over  the  equivalent  finite  diffe'^ence  model. 

Adding  an  extra  margin  to  the  already  established 
superiority  can  be  accomplished  by  using  the  low-resolution 
model  with  a  variable  grid.  Figure  44  is  similar  to  the 
tests  in  Figure  38  except  that  a  variable  resolution  grid 
(R  =  2.0)  was  utilized.   A  96-h  forecast  for  a  variable 
resolution  grid  (R  =  3.0)  is  shown  if;  Figure  45.   Comparison 
of  Figures  38,  40,  44,  and  45  show  tie  convergence  of  the 
low-resolution  uniform  solution  towards  the  high  resolution 
solution  as  the  degree  of  variability  is  increased.   Figure 
46  is  a  graph  of  geopotential  amplitudes  at  hour  96  versus  y 
mode  wave  number  (x  wave  number  one)  for  the  uniform, 
variable,  and  high  resolution  models.   The  R  =  2.0  case 
better  represents  the  control  than  the  uniform  low 
resolution  forecast  especially  at  low  wave  numbers. 
However,  at  wave  number  three  and  above  all  models  have 
similar  differences  from  the  control. 

The  harmonic  analysis  demonstrates  an  improvement  with 
increasing  variability  but  is  not  conclusive.   Another 
method  to  evaluate  the  improvement  with  increased  varia- 
bility is  to  look  at  the  difference  charts  between  the 
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Fig.    44.      As    in    Fig.    38    for    R    =    2.0 
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Fig.    45.      As    in    Fig.    44    for    R    =    3.0 
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control  and  each  level  of  variability.   The  difference  chart 
in  Figure  47  is  the  control  96-h  forecast  (high  resolution) 
minus  the  low-resolution  96-h  forecast.   Figure  47  shows  a 
maximum  difference  in  the  center  of  the  vortex.   Figure  48 
is  the  comparison  of  the  high  resolution  96-h  forecast  and 
the  R  =  2.0  variable  grid  96-h  forecast.   The  maximum 
difference  has  decreased  from  the  uniform  case.   Figure  49 
is  the  comparison  of  the  high  resolution  96-h  forecast  and 
the  R  =  3.0  variable  grid  96-h  forecast.   The  difference 
centers  at  the  right-hand  edge  in  Figures  47,  48  and  49  are 
a  result  of  post  processing  interpolation  errors.   Although 
the  maximum  difference  in  the  center  of  the  vortex  is 
greater  in  Figure  49  than  Figure  48,  the  latitudinal 
difference  gradient  is  less.   ,^oot  mean  square  differences 
from  the  high  resolution  control  are  652.2  m^/s^,  608.8 
m^/s^  and  469.7  m^/s^  for  the  low  resolution,  R  =  2.0  and 
R  =  3.0  tests  respectively.  The  forecasts  have  definitely 
been  improved  at  each  level  of  increased  resolution. 

G.    EXPERIMENT  4 

This  experiment  will  demonstrate  a  moving  grid 
capability.   Harrison  (1973)  demonstrated  the  use  of  a 
moving  grid  in  his  finite  difference  model.   The  inherent 
problems  with  moving  finite  difference  models  are   abrupt 
changes  in  resolution  causing  spurious  noise.   The  noise, 
unless  damped,  is  magnified  if  left  unattended  every  time  a 
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Fig.  47.   Forecasts  and  difference  charts  at  hour  96  for  the 
uniform  high  (S258115)  and  low  (S158115)  resolu- 
tion GFEM  models.   Perturbation  =  25.0  m/s. 
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Fig.  48. 


As  in  Fig.  47  for  the  high  resolution  (S258115) 
and  the  variable  (S178115)  (R  =  2.0)  GFEM  models. 
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Fig.  49.   As  in  Fig.  47  for  the  high  resolution  (S258115) 

and  the  variable  (S198115)  (R  =  3.0)  GFEM  models, 
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grid  moves.  In  a  GFEM  model  there  is  a  smooth  transition 
into  the  fine  grid  area  with  no  attendant  noise.   Therefore, 
moving  the  entire  grid  with  the  fine  grid  is  simply  a  matter 
of  accurate  interpolation. 

Three  tests  will  be  performed.   First,  uniform  low  and 
high  resolution  models  will  be  Integrated.   These  models 
will  not  use  a  moving  grid.  Then  a  model  with  the  same 
degrees  of  freedom  as  the  low-resoljtion  model,  but 
employing  a  moving,  variable  grid,  will  be  run.   The  high- 
resolution  model  will  act  as  the  control.   Identical  initial 
conditions  will  be  used  for  all  three  tests.   These  condi- 
tions are  similar  to  those  developed  for  experiment  3, 
except  that  the  source  term  will  be  zero  throughout  the 
course  of  the  integration.   The  ini;ial  conditions  include  a 
sharp  trough  with  a  relative  vcrticty  maximum  which  is 
advected  by  the  mean  flow  in  the  absence  of  any  forcing  (the 
source)  . 

Of  particular  concern  is  the  generation  of  noise 
incurred  during  grid  movement.   Therefore,  a  forecast  of 
96-h  is  made  during  which  this  grid  moved  eleven  times.   As 
previously  stated,  the  model  utilizes  absolutely  no  diffu- 
sion, filters  or  friction.   If  spurious  modes  are    excited 
during  the  movement  of  the  grid,  they  will  appear  in  the 
harmonic  analysis. 
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The  procedure  of  moving  the  finite  element  grid  is 
achieved  by  displacing  all  the  nodes  one  unit  distance  in 
the  direction  desired.   This  retains  the  relative 
distribution  of  the  nodes  and  merely  moves  the  entire  grid 
as  an  entity.   Ramifications  of  this  approach  are: 

1.  The  areas  of  each  element  remain  the  same  and  the 
numerical  quadrature  coefficients  for  all  the 
associated  integrals  do  not  have  to  be  recomputed. 

2.  The  unit  distance  moves  should  be  the  smallest 
grid  increment;  and 

3.  Direct  solvers  that  utilize  a  preprocessing 
procedure  based  on  the  grid  distribution  will  not 
require  the  preprocessing  to  be  recomputed  after 
each  move. 

The  most  important  aspect  of  moving  the  grid  is  to  accurate- 
ly interpolate  the  new  field  values  for  the  new  grid  point 
position.   In  this  dissertation  a  standard  International 
Mathematics  and  Statistical  Library's  (IMSL)  bi-cubic  spline 
interpolating  routine  was  employed.   The  movement  will  occur 
only  in  the  x  direction  and,  therefore,  a  one-dimensional 
interpolator  suffices.   The  interpolator  accommodates  the 
periodic  boundary  condition.   However,  interpolators  for 
non-periodic  boundary  conditions  are  available.   The  accu- 
racy of  the  interpolator  is  exact  at  grid  points  for  a 
uniform  distribution.   Initial  experimentation  verified  that 
a  field  could  be  moved  with  no  loss  in  accuracy  for  a 
uniform  grid.   Experiment  4  will  test  its  capability  to  move 
a  variable  grid  with  R  =  2.0. 
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H.         RESULTS    4 

Figures    50    and    51    are    the    initial    and    96-h    forecasts 
for   the   uniform    low-resolution   model.      The    perturbation    is 
25.0  m/s,  mean  flow   10  m/s  and  mean  depth  5  km.      The  initial 
divergence   is    a    zero    field   but    spins    up   through    geostrophic 
adjustment.      The    northeast/ southwest    trough    orientation    is    a 
manifestation   of  the   nonlinear    interactions    occurr- ng    for 
this    highly    active    perturbation.       Figure    52    is    the   96-h 
forecast    for   the    uniform    high-resolution    model.      The 
northeast/southwest    trough    orientation    is    evident    in   this 
forecast    as    well.      Figure   53    is   the   96-h    forecast    "or   the 
low-resolution,    moving    variable    grid    (R    =    2.0)    model. 

Time    steps    were   4320s,    3085.7s    and    3600s    for   ■;he    low, 
high    and    variable    models,    respectively.      The    variat)le    model 
took  95  full   time  steps  and  two  half  steps  during  1;he  96-h 
forecast.       The    half    steps    are    required    only    to    sta-t    the 
model   to  obtain  the  N  and  N  +  1  time  levels  in  t^e  leapfrog 
time    stepping.      The   grid    moved    11    times   during   the   forecast 
at   time   steps   8,    16,      24,    32,   40,   48,    56,    64,    72,    80    and    88. 

The   X    axis    origin    in    Figure    53    is   3530    km    indicating 
eastward   grid   movement.      The   vortex    is    still    in    the    fine 
mesh    area.      Additionally,    there    is    no    noise    in    the    forecast 
fields.      To    reiterate   a    previous    statement,    there   are   no 
filters,    diffusion    or   friction   terms   employed   during   this 
integration.      Comparisons    among    Figures    51,    52    and    53    show 
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Fig.    51.      As    in    Fig.    50    for    a   96-h    forecast 
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Fig.  52.   As  in  Fig.  51  for  a  high  resolution  GFEM  model 
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Fig.  51  due  to  the  movement  of  the  grid  during  the 
forecast. 
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that    the   variable    grid    forecast    more   closely   approaches   the 
control    than   does   the   uniform   low-resolution    model.      This    is 
to    be    expected    since   the    fine   mesh   area   has   been    moved    to 
coincide   with    the   act" ve   region    of   the    forecast   domain. 

Table   6    shows    th(2   geopotential    harmonic   analysis    for 
the   three   models    f.fter   a   96-h    forecast. 


TABLE    6 
Harmonic    Analysis;    (Geopotential)    for   Experiment   4 

Units    are   m^/s^ 


Wave 

Low  Re  sol 

1 

876.1 

2 

167.4 

3 

122.5 

4 

54.7 

5 

20.3 

6 

5.1 

904.7 

152.9 

176.9 

39.0 

10.2 

7.6 


Movi  ng    Grid 

850.0 

146.6 

159.6 

39.0 

9.6 

7.3 


Comparisons   of   the   higher   wave   number  amplitudes   show 
almost   identical    amplitudes   in  y   modes   four,    five    and    six 
for   the    moving    case    and    high-resolution    uniform    case.      The 
amplitude   in   y   mode    six   indicates   that   the   movement    of   the 
grid    has    not    generated    unwanted    noise.      The    overall    improve- 
ment   for   the    moving    case    is    a    result    of   the    finer   grid's 
ability    to    better    resolve   the   small-scale   interactions 
occurri  ng . 
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VI.   SUMMARY  AND  RECOMMENDATIONS 

Numerical  weather  prediction  requires  that  the 
equations  governing  atmospheric  flow  be  solved  numerically. 
These  flows  cover  a  large  range  of  scales.   Classical 
modeling  schemes  have  successfully  forecast  the  longer 
scales  and  improvements  are  being  sought  for  the  smaller 
seal es. 

In  this  dissertation,  the  GFEM  has  been  employed  for 
the  shallow-water  equations  in  a  channel  domain.   The 
differentiated  form  of  the  equations  on  an  unstaggered 
finite  element  grid  with  a  semi-implicit  time  scheme  has 
been  utilized.   An  equivalent  finite  difference  model 
representing  the  conventional  approach  has  also  been 
utilized.   Analytic  initial  conditions  representing  simple 
atmospheric  waves  and  a  more  complex  source  term  have  been 
used  to  test  the  different  models. 

The  hypothesis  of  this  dissertation  is  that  the  GFEM  is 
a  viable  option  for  numerical  weather  prediction.   The 
theoretical  foundation  of  the  GFEM  leads  to  a  minimization 
of  the  error  between  the  actual  equations  and  their  approxi- 
mation.  This  minimization  implies  a  basis  for  expected 
improvement  in  forecast  capabilities.   Three  features  were 
addressed  to  support  that  hypothesis.   First,  alternatives 
for  variable  grids  were  investigated.   Second,  forcing  near 
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the  smallest  grid  scale  determined  if  the  GFEM  model 
resolved  the  smaller  scales  better  than  a  finite  difference 
counterpart.   Finally,  a  proof  of  concept  for  the  ability  to 
move  the  grid  was  demonstrated.   Each  task  adds  to  the 
usability  of  the  GFEM  and  illuminates  the  well-rounded 
potential  that  it  contains. 

First,  the  ease  with  which  the  method  allows  variable 
grids  was  demonstrated.   In  all  cases,  no  numerical 
techniques  were  employed  to  damp  or  filter  any  field.   The 
grid  refinement  obtained  when  using  a  rectangular 
subdivision  was  more  attractive  than  with  a  triangular 
subdivision.   It  allowed  refinement  to  be  easily  obtained 
with  increased  accuracy  from  the  higher  order  polynomial  and 
decreased  computational  effort  due  to  reduced  operation 
counts.   Additionally,  it  was  shown  that  the  geometry  of  the 
triangular  elements  greatly  influences  the  problem. 
Furthermore  unless  specific  constraints  are  enforced,  the 
boundary  conditions  would  not  be  satisfied. 

Second,  the  ability  of  the  GFEM  to  resolve  atmospheric 
phenomena  occurring  near  the  scale  of  the  smallest  grid- 
length  was  shown.   A  small-scale  steady  state  solution  was 
cast  into  a  source  required  for  that  solution.   Integrations 
showed  that  the  GFEM  model  resolved  the  small-scale  forcing 
admirably  and  far  exceeded  the  performance  of  an  equivalent 
finite  difference  model.   Even  when  the  resolution  of  the 
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finite  difference  model  was  doubled  and  tripled,  it  could 
not  match  the  performance  of  the  low-resolution  variable 
grid  GFEM  model.   The  responsiveness  of  the  GFEM  model  near 
the  smallest  resolvable  grid  length  clearly  demonstrates 
that  the  variable  grid  concept  is  worth  the  computational 
effort. 

Finally,  the  concept  of  moving  the  variable  grid  was 
demonstrated.  The  ability  to  resDlve  a  small-scale 
phenomenon  is  a  tremendous  asset  to  the  GFEM.   The  ability 
to  move  the  grid  with  the  atmospheric  phenomena  further 
enhances  its  usefulness. 

Throughout  this  dissertation  the  GFEM  model  has 
demonstrated  desirable  characteristics.   It  has  conservative 
properties.   It  propagates  atmospheric  waves  better  than  an 
equivalent  finite  difference  model.   It  allows  variable- 
resolution  grids  and  responds  better  than  an  equivalent 
finite  difference  model  near  the  smallest  gridlength. 
Moving  grids  can  be  achieved  with  no  apparent  noise 
generation.   It  can  utilize  direct  solvers  and  is  a  natural 
choice  for  vectorization  on  large  computers.   It  truly  is  a 
viable  option  for  simulation  of  atmospheric  flows. 

The  success  of  this  dissertation  mandates  further 
research.   The  next  logical  step  would  be  a  baroclinic 
model.   Completion  of  a  baroclinic  model  should  allow  its 
application  as  a  regional  model.   Additionally,  different 
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basis    functions   could    be   employed    for   the   vertical 
discretization,    such    as    solutions    to   the   vertical    structure 
equation.      For   the   moving    variable-resolution    grid,    further 
research   into   two-dimensional    movement   is   required.      The 
proof   of   concept    in    this   dissertation   of   one-dimensional 
movement   can   logically   be   extended   to   two   dimensions.      Upon 
completion   of  two-dimensional    movement,    an   operational 
forecast    capability   for   tropical    cyclones    should    be    investi' 
gated.      The    superior    small-scale    response    of   the   GFEM 
indicates    potential    increase    in    skill    for   both    regional    and 
tropical    cyclone/hurricane    prediction. 
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APPENDIX  A 
STABILITY  ANALYSIS 

A   one-dimensional    linear   analysis    of  the   forecast 
equations    will    be    presented.      The   analysis    was    obtained 
through    personal    communication   with   Dr.    Andrew   Staniforth 
(Recherche   en    Prevision    Numerique).      The    primitive    form    of 
the   forecast   equations    and   a    semi-implicit   time   scheme   will 
be   analyzed,    because   the    results    are    identical    for   the 
vorti ci ty/di vergence    form.      The   one-dimensional    equations 
wi  th   a  mean   f  1  ow   U    are 

at     3x  ax  i^\-i; 


at  ax.  ^'^  ^  I 


li   +  c^iy.  =  -  u  ii  (A-3) 

at       ^  3x  \i^  ^) 

Time   derivatives    are   evaluated   with   a   centered    time   dif- 
ferencing   and    the    other    terms    on    the    left-hand    side   are 
averaged    between    time    levels    (t    +    At)    and    (t    -    At). 
Equations    (A-1),     (A-2)    and    (A-3)    become 

U(x,t+At)    -   u(x,t-At)         1  r((j)(x+AX,t-l-At)    -    (t>(x-AX,t+At) 

2At  2  L  2ax 

+    ('l)(x+AX.t-At)    -   4)(x-AX,t-At)-|    _        r.    rU(x+AX,t)    -    u(x-Ax,t)-, 
2ax  J  -  -  u  l  2ax  -• 

+  f  v(x,t)  (A-4) 
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v(x,t+At)   -  v(x,t-At)  r,  rV(x+AX,t)   -  v(x-Ax,t)-[ 

2At  -  -  u  L  2ax  -■ 


-     ^    'j(X.t)  (;^.5) 


(<)(x,t+At)    -    (f)(x,t-At)     .    $    r(uCx+AX,t+At)    -    uCx-AX,t+At) 

2At  2  ^  2  AX 


^    Cu(x+AX,t-At)    -    u(x-AX,t-At))n  H  /  4*  C  <+AX ,  t )    -    ({)(x-AX  ,  t)  x  /a     ft\ 

2ax  J  ~   -  ^K  2ax  j  vm    o; 


Assuming    that    a    function,    F,    can    be    expressed    as 

F(x,t)    =    F'e^^*^^^   ■'w^)  (A-7) 

then 


F(x,t..t)  ^  F(x.t-At)   .  ,.   ,,3(,,,,;<^'<  ^  «*'  (A.9) 


Equations  (A-4),  (A-5),  and  (A-6)'  become,  after 
substi  tuti  on , 

1^'  -  fv'  +  ikU'  =  0  (A-10) 

At 

fu'  +  ^'   =  0  (A-11) 

ikUu'    +^'  =  0  (A-12) 
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k  A  X 

where  s  =  sin  wAt  +  k'uAt,  c  =  cos  wAt  and  k'  =  sin  —rz-» 
The  determinant  of  the  system  of  equations  (A-10),  (A-11) 
and  (A-12)  when  set  equal  to  zero  yields 


At 


(A-13) 


The  roots  are  s=0  and 


s^  =  CfAt)^  +  c^UAt)^  $ 


(A-14) 


The  roots  of  this  second  order  equation  when  requiring  w  to 
be  real  yield  the  stability  criterion 


At  < 


1 


JJ 

AX 


+  f 


(A-15) 
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APPENDIX  B 
NUMERICAL  QUADRATURE 

Evaluation  of  the  Galerkin  integrals  requires  a  fast, 
efficient  method  .   Two  different  elements,  triangles  and 
rectangles,  have  been  employed  and  the  quadrature  schemes 
for  evaluation  are  presented  here.   In  the  first  case, 
triangles,  area  coordinates  are  used.   For  the  rectangles, 
integration  formulas  are  basF.d  on  an  orthogonal  axis  trans 
formation.   In  both  cases,  the  integrals  to  be  evaluated 
contain  either  products  of  basis  functions,  products  of 
derivatives  of  basis  functions,  or  a  mixture  of  both. 

Zienkiewicz  (1971)  developed  the  relationship  between 
the  Cartesian  coordinates  of  a  triangle  and  the  area 
coordinates  as 


X  =  L^X^  .  4^2  '  hh 


(B-1) 


Y  =  L,Y,  +  L^Y,  +  L.Y 


•ri 


3"3 


(B-2) 


i  =  h  ^  4  ^  4 


(B-3) 


where    (x2,yi),     (X2,y2)    an<l    (^s.ys)    a>"e    the    Cartesian    coordinates 
of    the    triangle's    vertices    and    L^^,    L2    and    L3    are    the    area 
coordinates.      Here 


L^   =  A^/A 


L2  =  A2/A 
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4  =  A3/A 


where  Aj^,    A2    and  A^    are  shown  in  Figure  54  and  A  is  the  area  of 
the   entire   triangle.      He    further    shows   that 


// 


I    3  I    b  ,    c    ,    .  alblcl         on 

h     4     4    '"'"y  -  (a^b+cf2):  2A 


The  L's  are  the  basis  and  test  functions  used  during  the 
approximating  process.  Figure  55  defines  some  necessary 
parameters    for    evaluation    of   the   derivatives. 

Equations  B-1,  B-2  and  B-3  can  be  written  in  matrix 
form   as 


hi 


L^3J 


J_ 
2A 


2A 
2A 
2A       b 


'1 


^2       ^2 


'3J 


ri 
X 
Y 


(B-5) 


Differentiation  of  Eq.  (B-5)  shows  that 


3  b. 


3X 


i=l  2A  3L. 


(B-6) 


^-     2 


3  a 


i  3 


ay 


.  T  2A  dl. 

1=1    1 


(B-7) 


However,  the  derivative  of  the  basis  function  with  respect 
to  the  area  coordinate  is  non-zero  only  where  the  basis 
function  and  area  coordinate  coincide.   In  fact,  the  non- 
zero value  is  one.   Therefore,  using  V-j^  as  a  basis  function 


9X    2A 


CB-8) 
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Fig.  54 
Y 


Cartesian  coordinates  vs.  natural  coordinates 


i 


h    =  ^3-  ^1 


.83   Xj-  x^ 


b  =  y  -  ¥ 

1    '2         ': 


2     '3 


^  =  Yi  -  Y: 


-*-X 


Fig.  55.   Triangle  definitions  for  area  coordinates. 
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All  integrals  can  be  evaluated  by  combination  of  these 
quadrature  formulas  for  triangles. 

An  orthogonal  axis  transformation  will  allow  similar 
quadrature  formulas  for  rectangles.   The  rectangle  shown  in 
Figure  56  can  be  transformed  by 


x-x. 


y-y. 


c  = 


n  = 


(B-9) 


The  values  of  ;  and  n  at  each  corner  are  shown  in 
parentheses.   The  basis  function,  V^  ,  becomes 

Cl+?,.d(l+n,n) 


V.  = 
1 


(B-10) 


Fig.  56.   Orthogonal  axis  transformation  for  rectangular 
integration  formulas. 

Derivatives  of  the  basis  function  become 

3V.    1  3V. 


3X 


a  3c 


and 


3y  ~  b  3n 


(B-11) 
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By  substituting  and  integrating,  the  4  x  4  matrix  with  the 
interaction  coefficients  can  be  determined  for  a  derivative 
or  straight  inner  product.   For  instance,  the  straight  inner 
product  i  s 

1  1 

dcdn 


:.  .  =//v.V.dxdy  =  ab  ^  ^  V.V. 


1  1 


=  T|(4s-^jH24n,n,) 


(B-12) 


and   the  mixed   derivative   is 

aV.   aV 


''      ■   3 v..   9V 


-^l   ^i^j^^/^    (Hn.n:.(Hn.n)dn 


l|^C2c,Cj)(24n,nj) 


(B-13) 


These  quadrature  rules  allow  the  evaluation  of  any  integrals 
when  using  rectangles. 
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